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FOREWORD 


Stanford University and the University of California at Los Angeles 
record in this Final Report, the results of a nine-months study on the 
application of a Drag-Free Satellite to geodesy. The work was performed 
under NASA Grant NAS-12-695 which was awarded in response to our proposal 
[Refs. 1 and 2] "To Prepare a Preliminary Design of a Drag-Free Satellite 
and Study Its Application To Geodesy," and the Addendum dated January, 
1968. The work reported is principally in two areas: the first is the 

feasibility of making geophysical measurements which are not possible 
with conventional satellites, and the second area, preliminary design 
work on attitude and translation- control systems for spinning vehicle 
and possible coupling of attitude and translation control for gravity 
stabilized vehicles. 


A Drag-Free Satellite control reference is an unsupported proof mass 
shielded from all external surface forces by the satellite. Since only 
gravitational forces act on the proof mass, it follows a purely gravita- 
tional orbit. A control system in a satellite senses the relative posi- 
tion of the satellite with respect to the proof mass and actuates reaction 
jets forcing the satellite to follow the proof mass. The satellite there- 
fore also follows a purely gravitational orbit [Refs. 3 and 4], This 
concept has been developed to a high degree at Stanford under NASA Grant 
NsG 582 and was independently proposed at UCLA for aeronomy studies in 
1962 [Ref. 5]. There are two principal advantages of a Drag-Free Satel- 


lite for geodesy. 


A Xi O i, m 


sll i+r ~ 


cancels surface forces which might 


mask the small variations in the gravitational field; hence, it protects 
the long period perturbations due to weak effects from distortion and 
allows them to become large enough to be measured. Secondly, the perigee 
of the satellite can be made lower without introducing additional distur- 
bances to the satellite epheineris. It is possible, therefore, to enhance 
the effect of higher gravitational harmonics. 


Two types of geodesy information have been considered in detail. 

The first is an improved analytical model of the effect of tidal forces on 
the earth. Latitude and longitude dependence of the elastic deformation 
and energy dissipation have been included. From this mathematical model, 
it is hoped to determine accurately, the tidal interaction of the earth, 
moon, and sun. Some fixed tesseral harmonics of the earth's gravitational 
field have been measured by conventional satellites. A large number of 
tesseral harmonics (72 pairs) were selected as being uncertain and of 
interest. Through careful analysis, four orbits have been selected in 
which a Drag-Free Satellite can identify 69 of these 72 pairs of tesseral 
harmonics. The necessary theory for this evaluation has been reviewed and 
the analysis leading to the selections of these four orbits is presented 
in detail. 


In the design of the Drag-Free Satellite, there are two different 
types of performances that concern the designer. The control system 
which thrusts the satellite so that it follows the proof mass must be 
efficient in the utilization of propellant. Low propellant expenditure 
and tight control are the most important, measures of control performance. 
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A very different type of performance is discussed in evaluating how closely 
one achieves a purely gravitational orbit. The proof mass may be perturbed 
by a sensor, gradients in local magnetic field, mass attraction from the 
satellite, and many other phenomena. Mass attraction of the satellite 
on the proof mass is the largest of these very small forces which prevent 
the orbit from being purely gravitational. When the average value of these 
disturbances are along- track, they produce significant perturbations. 

To reduce this effect by several orders of magnitude, the satellite may 
be spun. Therefore, a thorough preliminary design study was made on a 
spinning vehicle, attitude and translation control system for a drag-free 
operation to establish the feasibility of achieving this significant im- 
provement of performance in case it is needed. The importance of adding 
a Drag-Free Satellite control system to an existing satellite was also 
recognized and hence, feasibility studies were performed to insure that 
special coupling effects would not degrade the translation control system 
performance. These are presented herewith in detail. 
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I 


TIDAL FRICTION WITH LATITUDE- DEPENDENT AMPLITUDE AND PHASE ANGLE 


Tidal disturbing functions were developed in which the ampli- 
tude factor k and lag angle c are expressed as sums of zonal 
spherical harmonics* 

In regard to the current evolution of the moon's orbit, the 
existence of a second degree harmonic in the lag angle could make 
a significant contribution to energy transfer to the moon; it is 
unlikely, however, that it has an important effect on the overall 
time-scale of the orbit evolution. 

If the moon formed in an equatorial orbit about the earth, 
tidal friction could do nothing to incline the orbit. Once the 
orbit was inclined, tidal friction could increase the inclination 
further if there was a commensurabi 1 i ty between the earth's 
rotation and the moon's revolution. 

Latitudinal variations in tidal properties would have 
appreciable effects on close satellite orbits- An appreciable 
second degree harmonic in the phase lag is needed, however, to 
reconcile the available data with the rate of the earth's 
rotational decel erat ion . Further progress requires satellites 
free from surface force effects* 
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A. INTRODUCTION 


As has been established by various studies in recent years, 
integration of the moon's orbit back in time under the influence 
of tidal friction with the dissipation factor 1 / Q inferred from 
observations of lunar 8 solar motion bring the moon close to 
the earth only 1.0 to 2.0 x 10 9 years ago (Gerstenkorn , 1955; 
MacDonald, 1964; Goldreich, 1966). Furthermore, at the closest 
point the moon's orbit has an appreciable inclination with res- 
pect to the earth's equator, which would rule out a fission 
origin. Goldreich (1966) and Munk (1968) suggested, however, 
than an asymmetric tidal response of the earth may have had 
a significant effect on the orbital evolution of the moon. 

Additional data which have recently been developed use 
the effects of tidal potentials on artificial satellite orbits* 
Kozai (1968) analyzed perturbations of the inclination of three 
satellites of 33° to 50° inclination. He obtained Love numbers 
k 2 varying from 0.23 to 0.33 with uncertainties less than + .03, 
and phase lags from 0° to '9° with uncertainties of + 5 to 7° * 
Newton (1968) analyzed perturbations of the node and inclination 
of four polar satellites- He obtained Love numbers k s varying 
from 0.28 to 0.44 with uncertaini tes less than + .03, and phase 
lags from 0*7° to 2-5°, with uncertainties less than + 0.8°. 

Newton emphasized that there appeared to be significant rate- 
dependence of both amplitude and phase angle- However, we still 
might hope to explain some of the differences among artificial 
satellite analyses, as well as the discrepancies from the 2° - 2-5° 
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lags inferred from the deceleration of the earth's rotation 
by an Asymmetric tidal response: i-e., one which is a function 

of location within the earth. 

B. DEVELOPMENT OF DISTURBING FUNCTION 


Let the disturbing function of the sun or moon be repre- 
sented in terms of the Kepler elements of the disturbing body 
and the spherical coordinates of the point of calculation of 
the potential (Kaula, 1964): 

U = 1 B tn C Smpq ^ P 'ta (5in ' [,) 

£,m,p,q 

• 

even f V 

* i c ? s } {v* - m(X+0)f, 

odd 1 ^ J 

(1) 

where r, cp, X are radius, latitude, and longitude 
9 is the Greenwich sidereal time; P^Csincp) is the 
Associated function; and 

respect i vel y ; 
Legendre 

B tn, ” Gm 4rasK( 2 ' 6 °"') 

(2) 

C tmpq - W*> W**' 

(3) 

= (4-2p)w* + (£-2p+q)M* + rnfl* 

(4) 
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where G Is the gravitational constant; m* , a * , e* , I*, M* , <ju* , 
ft* are the mass and Kepler elements of the disturbing body; 
and p£ m p(l*) and (e*) are polynomials as derived by Kaula 
( 1966 ) . 

Now, regardless of the nature of the tidal response of 
the earth -- oceanic, bodily, or otherwise -- so long as there 
is no significant non-linearity (i*e., tidal terms interacting 
with themselves or other terms), then the tidal potential at 
the earth's surface r = R can be written as: 


T(R. <P, X) - l k t (<p. X) R l B £ m 

-tmpq 


pq 




-t-m even 
t-m odd 


{ v tmpq ' € ^npq (tp ’ X) " m ( x+e )} < 5 > 


We should expect the Love number (cp , X) and phase lag 
e^mpqC^* M to be rather smoothly varying functions of position, 
and hence representable by spherical harmonics; this would cer- 
tainly be true of that portion of their effect which would affect 
satellite orbits* The full development of Eq. (5) would Involve 
products of tesseral harmonics, and hence would be most con- 
veniently done in complex representation, along the lines of 
Kaula (1967)- But a fairly elementary consideration shows that 
any longitude-dependent part of the product will have 

its effect on any orbit "averaged out" by the earth's rotation. 
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Hence for the purpose of orbital analysis let us assume and 
e^pq to be functions of latitude only, which allows us to retain 
a real representation: 




y H iu p ho (sjncp) 

h 


( 6 ) 


and 


e< tmpq = y C n(tmpq) P no^ s ' nC! ^ 
n 


( 7 ) 


In derivations where the frequency dependence of e n is not being 

emphasized, we shall drop the argument (-tmpq) • 

If it is assumed that all the e are small enough that 

n 

sIn ( £ -tfflpq) = € > cos ( C <mpq) = ' • then the tidal P otentia ' at 
the surface can be written: 

T < R > I K tmpq K th p h 0 <siTO) P t™ (slncp) 

-tmpqhn 



l-m even 
-t-m odd 




. -s l-m even 
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where 


X 


^ ^mpq 


(9) 


The product P^, m (sincp) P irt (sincp) can be converted to a sum: 


jo 


Uj 

^jkm- ^km 


( 10 ) 


Values of m are given in Table 1* Applying Eq* (10) to 
all products in Eq* (8) allows us to write the tidal potential 
T(r) at any radius r > R as: 


T(r) " I K 'tmpq K th Q ^hkm 
-tmpqhnks 


• [(?)“' ’Jsaz r (’w, - •»*«) 


+ 6 o /jR") p 

n y knsm\r/ r sm 1-cosJ l v 2mpq 


s+1 


f.-L 'I't-m even r 
,{- S cos}. . . . K 


m ( \+ 6 ) 
(11) 


>] 


Using the usual conversion from spherical coordinates to Kepler 
elements of the perturbed satellite orbit (Kaula, 1966, p. 37), 




f cos 
s 1 n 
even 


t 


i 

i 

( -t-m) 


s 1 nl 
-cosj 
odd 


(T2) 
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we get 


^ Z ^-tmpq H ^h ^<hkm 

-tmpqhnksj g 


k+1 me (k-^)even 

W> W>{f.,A,n). 


(k-£) odd 


•{ 


v kmj g ^mpq 


} + E n <J 

( s-l ) even 


s+1 


knsrrAa^ ^smj ^ 1 ^ 


-sin x 7 vtvu 

G sjg^ e) 4" 1) m cos} (s-<t) Qdd { v smjg " 'tmpql 


03 ) 


To simplify, let us substitute (k^e)^ for Hh Z_ € n Q 
then we can write 


n 


knsm' 


k+1 


T = 


Z ^-tmpq(a) ^kmj ^ ^kjg^ Q-thkm 

-tmpqhkj g 


• [»th{(-l7sin} (k . t)odd + ( k ^)h{(-,rcos} 


cos . even / >, -sin (k-*) even- 


(-1) cosJ (k . t) odd 


{a - v* } 
l kmjg -tippqJ 


0*0 


For analysis of the moon's orbital evolution, we require 
(1) zero rate: 


^kmjg '^tmpq 


05 ) 
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or , from Eq • ( 4) , 


j = £(k--t) + p 

• 

( 16 ) 

g = q; 

and (2) small "damping" factor ^ 

la* 4+1 a k+, )' , ( or 

(17) 

Ik = 20, 22, 24, 

33, or 42 

( 18 ) 


The combination 31 for Ik is excluded because a first degree 
harmonic would represent a shift of the origin, the earth's 
center of mass . For lunar orbit evolution then let us write 


T - T, r ' 


0 T ! 2 ' r *4 


where 


T 0 “ I K 201q(a) F 000 G 00q Q 2h00 *201 hq 




qh 


^2 ^ ^2mpq(a) F 2mp G 2pq ^2h2m *22mphq 

qhmp 


^4 = ^ I ^-tmpq(a) F (6--f-)mj G 2Jq 

1=2 qhmp 


(19) 


> 


^( 6-t) h-tm * (6- t)^mphq 


J 
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in which 


♦tkmphq = ( K th cos * ( k '. e )h si ") ( v km j q ‘ 'tmpq) 

and j is related to k, l, and p by Eq. ( 1 6) . 

For analysis of tidal perturbations of artificial satel- 
lites, we require (1) no short period terms, or 


9 ° 2j - k (20) 

and (2) small "damping" factor , or 

l - 2 ' ( 21 ) 

Eq. ( 1 4) thus becomes 


V /'R'\ k+ ^ 

L . ^2mpq\a/ ^kmj ^kj(2j-k) ^2hkm 
mpqhkj 



k 

/ C0S \ 

even 

j -sin 

even" 

• 

*2hl(-1) m sini k 

+ \ k t e )h 

odd 

l(-1) m cos/ k 

odd .. 


- [ V kmj(2j-k) ‘ 'ImpqJ < 22 > 



Values of the amplitude factors K 
are given in Table H- 


■4mpq 


and associated rates 
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c. EFFECTS ON LUNAR ORBIT EVOLUTION 


Of the terms in Eq./(19), those with ^mpq = 2010 are 
meaningless because ^gOOO ' s zero. The lowest value of h which 
gives a non-zero Q 2 h00 ' s 2 * frrom Table t“2, K 20 jj is appreciably 

larger than K 201-1* ^001 * s coe ^^' c ' ent °f cos M in the 
expansion of a/ r , or e. Writing out the parts of K 2 qjj and 

evaluating Q 2 200 ^ rom Table I we get 

T ° “ R3C^ slna| * • i] [| e * + 7S e * 3 + • • •] 



Hg s cos 


v.ao e 2 sin) (M- M* ) (23) 


From the Lagrangian planetary equat ions-of-mot ion , we get 
(dropping asterisks after differentiating To): 


a 0 


JL Ms. 
na SM 


Gm* 2e^RN 4 fi 
na s 5R'a/ lh 


s i n 


*] [| e + U e * + 





3Gm*e 

10na 2 Rj 


(I) 


Hao 


(24) 


and 
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eo 


1-e ^Tn 
na" J e ?M 


H^(|) [| s in a l - i] [Je + y£e 3 + . 


( l<8 e a ) 


3Gm*e /RV 
“ 20n^R\i/ * 20 €a 


(25) 


These terms express the effect of the interaction of the 
eccentricity of the moon's orbit with a variation of the tidal 
dissipation in the earth which is symmetric about the equator. 
Since the earth's rotation is not involved, there is no angular 
momentum transfer: only energy transfer. Eq- (24) states that 

if dissipation is concentrated near the poles in the earth, 
the moon will move away faster; while if it is concentrated 
near the equator it will move away more slowly. Eq. (25) is 
merely consistent with conservation of angular momentum, or 
8 * 0 - 6 *)*. 

Putting in numbers, we have for the present lunar orbit, 
from Eq. (24), in planetary units 

a 0 = 1.08 x 10‘ 13 *20 62(2011) (26) 

From Eq. ( 4l ) of Kaula (1964), 
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as ~-^(f) 6 Hao e 0 ( 2200) 

= 0.97 X 1 0"” 1 3 Kgo ®o ( 2200) (27) 

The frequency (in an earth reference frame) for £mpq = 2200 
is about 55 times that for -tmpq = 2011. Hence, since a sig- 
nificant part of the tidal dissipation is in the oceans, probably 

G 0 ( 201 1 ) « e 0 ( 2200) (28) 

Furthermore, we would expect e 0 + e a P 3 to be greater than zero 
everywhere, so probably 

~ c 0 ( 20 11) < € a ( 201 1 ) < 2G 0 (201 1 ) (23) 

Thus it seems unlikely that the effect on the lunar orbit 
evolution of non-uniformity of dissipation can explain the dis- 
crepancy between the lag angle Inferred from lunar & solar motions 
(most recently by Curott, 1966) and that from artificial satel- 
lite orbits by Newton (1968). Since is positive, for non- 

uniform dissipation to explain some of the discrepancy its 
contribution to should be negative and hence a 0 should be 
positive, corresponding to a predominance around the poles 
(see Kaula, 1968, p. 197) • 

The most significant result of Goldreich (1966) was that if 
dissipation has been uniform over the earth's surface, the 
inclination to the equator could never have been less than 10° 
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within 10 earth radii. Goldreich suggested that strong local 
dissipation in a few places in the earth's oceans or crust may 
lead to unanticipated deviations from his results* No matter 
how "local" was this dissipation, and no matter how close the 
moon, the main effect on the orbit must be expressible in spheri- 
cal harmonics -- just as is the case for satellite orbit per- 
turbations by gravity anomalies fixed in the earth. A faster 
rate of rotation such as prevailed early in the earth's history 
would make this all the more true. 

For the inclination to be changed by the disturbing function 
of Eq* (19) » the combination F tmp F kmj[ (k - 2 j )cot 1 - m csc '] 
must be non-zero (see, e.g., Kaula, 1966, p. 40). For I = 0, 
this never occurs: contains a non-zero factor only 

for m = k-2j . If m ^ k-2j , F^ mp F, <m j is at least order I 3 . All of 
which goes to say that purely latitudinal variations in tidal 
properties can do nothing to wrench the moon out of an equatorial 
orbit. What is necessary is an interaction such that 


[(k-2j) - (£- 2p) J (ui+M) + (s-m) (ft-8) 0 


(30) 


with s 3 ^ m: i.e*, a resonance. To maintain such a resonance 

• • 

long enough to have significant effect, ri/9 must have about 
equaled n/9, where n is the mean motion, related to u = G(M+m*) 
by Kepler's third law, 


n 



(31) 
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From conservation of angular momentum (Eq. (59) j Kaula, 1964): 

ce = Jjm * |i 8 / 3 .n“ 4 ^. n (32) 

where C is the earth's moment of inertia, or 

n/8 = 3cm*" 1 u" a/3 n 4/3 
- 3cm*" 1 a“ a 

« 81/a 8 (33) 


for a In earth radii. Integrating the moon back close to the 
earth gives a rotation period of about five hours (Fig. 9» 
Goldreich. 1966), or 0-28 for 9 in the "planetary" units of 
radians/807 sec- Setting n as equal to 0.28 x 8l/a 2 and equating 
it to n by Kepler's law (Eq. (31)), we get an absurdly large 
semi-major axis a. Hence it seems unlikely that any longi- 
tudinal variation in tidal properties ever had a significant 
effect on the lunar orbit. This holds true even for non- 


linearities in the tides, since in Eq* (30) the tidal harmonic 
ks now has a completely independent specification from -Cm. 

The next, more desperate, possibility is that an irregu- 
larity in tidal properties would affect n so much in a direction 
counter to the central tidal term as to hold the moon at a 


resonance point where its inclination could be changed. A 
situation might have existed which is mathemat ical 1 y similar 
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to that of Mimas and Tethys , with a tidal bulge taking the place 
of the inner satellite Mimas- However, the treatment of Allan 
(1968) indicates that the inclinations must have been non-zero 
before these satellites became coupled, and that the subse- 
quent increase in the inclinations depends on a factor of 
order \ s . 


We thus seem to be forced back to the conclusion that at 
least part of the moon was captured. The interesting question 
then is how small a portion of the moon needs to be captured, 
taking into account that resonance with longitudinal variations 
in tidal properties may help to further feed the inclination. 
Since the mixed capture-fission hypothesis would involve seve- 
ral more ramifications that either theory alone (see the dis- 
cussion & conclusions of Goldreich, 1966), it seems appropriate 

*- « ^ f i /-%*•> r»pr>o»* . 

tw MC I & l I M I l I »ci I V*wi O lUCI UHWII VMWI w V W.MN/ .WI " 

In regard to the time-scale problem of the lunar orbit 
evolution, the conclusion of Eqs* (26) - (27) that latitudinal 
variation in tidal lag can account for only a minor part of 
the current evolution indicates that it was of even less sig- 
nificance in the past, and hence a less important effect than 
changes in the extent of shallow seas, as suggested by the most 
recent paleontological work (Panel la et al • , 1968). 


D. EFFECTS ON CLOSE SATELLITE ORBITS 

We wish to examine the bahavior of artificial satellites 
of small a/R under the influence of the disturbing function 
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given by Eq. (22) with a view to explaining results already 
obtained and suggesting specification of future orbits to 
determine tidal properties of the earth. 

In his analysis, Kozai (1968) determined k.and e from the 
perturbation of the inclination A| of argument 0. Hence he 
used only the term mpqhkj = 110021 of Eq. (22), or 

T 2 ■ K 2tOo(f)’ '[' 2 sln 1 cos 0 0 -® a )' 3/S 

* K3o(cos -e 0 sin) (0-0*) (34) 

Newton (1968) determined k and e from the perturbations 
of inclination A| and node A0 with arguments containing 0 and 
20. He obtained the lunar and solar-orbit dependent factors 
by numerical Integration, which would be equivalent to using 
all 6 terms m = 1 ,2 and p a 0,1,2, with qhkj = 0021 in Eq. (22): 

T 2 " l K 2mpo(!) a <'-e 3 >- 0/S 

mp 

• Hs o ( cos -c 0 sin) m( 0- 0*) (35) 

Since the perigee argument iu is absent from the disturbing 
functions (34) and (35)* odd zonal variations h = 1,3,5 ••• • 

in k 3 and e could not have given rise to perturbations of the 

% 

same period as equations (34) and (35); only even variations 
h «= 2,4 . . . Hence in the expressions for perturbation of 
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A| (dependent on bR/ 20) containing the Love number k we can 
make the substitution 



*210 


k,c 




I (a) F kmk/ 2 G k(k/2)0 
k , h 



(36) 


m = 1,2; h,k even 


and in the expressions for the perturabtion AO (dependent on 
*R/8l), the substitution 



bF 


2ml, 


b I 210 



V ^ k+1 5fr kmk/2 P 

L W 81 b k(k/2)0 

k ,h 


Let us define 


^2hkm 


( . ^ 
j k n 

J( kg £ ) 


h 


i 


m = 1,2; h,k even 


(37) 


^hm ~ Zj \a/ F kmk/ 2 G k(k/2)0 ^2hkm 
k 


(38) 


and 



( 39 ) 
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Then (36) and (37) become 



Kozai (1968) used satellites of essentially two specifi- 
cations : a/R = 1*30, I “ 33° * e ** 0.1 6 ; and $/r = 1*22 + .04, 

I a 48.6° + 1.5°» e = .01. Newton (1968) used satellites of 
one specification a/R = 1*157 + *016, I = 90.2° + 0*3° , e < .008. 
Values of corresponding to the m values of terms used in 
these analyses are given in Table 1-3. 

In Table 1-4 are summarized the results obtained by Kozai 
and Newton. Newtonls values for k and k6 are based mainly on solar 
perturbations and are combinations of values ranging from 0*31 
to 0.38 which are weighted inversely proportionate to their 
radiation pressure perturbations* In combining Kozai's values 
for his 47*2° and 50.1° inclination satellites we have utilized 
a similar weighting: 1/-21 for the 47-2° (ECHO 1 ROCKET CASE) 

and 1/.07 for the 50*1° (ANNA IB). In defining the lag angle 
we have accepted the Darwinian assumption, used by Newton, that 
it is proportionate to frequency: i.e*, e «* 6 for m ® 1 and 

c = 2& for m = 2 • 

We also have included in Table IV the lag inferred from the 
earth's deceleration, using the rate derived by Curott (1966). 

It appears in the column kj6, since it depends on *R/ *M and 
the tmpq => 2200 term is dominant (Kaula, 1968, pp . 198-203). 
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Applying Eqs • (40) - (41) to the close satellite data in 
Table IV we get as observation equations: 

’ *310, ~ .086 , + . 042^ fxzo *(*36)0 
■ .429, - -204, + .126 
• 986, .155, -.304 

.986, - .008, - .026 J 




(42 


Which yield 


ks = .372 - -240 P ao * -013 P 40 


(43) 


(k6) 3 = .006 - .093 P30 ” -079 P40 


(44) 


The large *32 of - .240 makes significant terms » when 

both h and n are non-zero, appearing in Eq. (13). If we let 


l 1 

t,n 


tnko K H 6 2n " 


2k 


(45) 


we get 


6 a = .003 - .028 P 20 - .263 P*o, 


which is even more implausible- 


(46) 
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An alternative procedure would be to use Newton's data 
alone for nso and H 22 * and then use Newton's data with the moon 
for (k6) a0 , (k6) 22 , (kS)^ . We get 

k 3 « . 351 - -055 P 20 (47) 

(k6) a - .0072 - .0121 P 2 o - -0004 P 4Q (48) 

and 

*>a = .0195 - .0331- Pao - .0028 P 40 (49) 

Eq. (47) is much more reasonable than Eq. (43), and Eq- (49) is 
somewhat more reasonable than Eq. (46). Physical necessity requires 
only that 630 is positive; the existence of ocean tides makes it 
possible that one of the zonal coefficients can be larger than 6 30 . 
Eq. (49) says that beyond the latitude where Pao is .0195/ *0331, 58°, 
phase leads should predominate over phase lags in the tide meter 
readings. A statistical analysis would be worthwhile- A major 
effect by ocean tides would invalidate the Darwinian assumption and 
make it impossible to extrapolate any inferences to rates other than 
semi-dirunal S- diurnal, such as the monthly rate which Eq. (26) sug- 
gests might be important in lunar orbit evolution- 

With the data on hand we get only a tentative indication 
of significant latitudinal variation in the earth's tidal pro- 
perties- But it is certainly a good enough indication to war- 
rant further effort in two directions: (1) the placing in 
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orbi t of artificial satel 1 i tes which are i ns t rumen ted to com- 
pensate for any surface forces, in order to determine long- 
periodic (fortnightly or more) variations in the orbit; and 
(2) the transformation of tidal models of the earth to a form 
compatible with the potential of Eq* (13) or ( 1 4) . 
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TABLE 1-1. FACTORS FOR PRODUCT TO SUM CONVERSION 


OF LEGENDRE ASSOCIATED POLYNOMIALS* 
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TABLE 1-3 . l&TLUENCE COEFFICIENTS OF TIDAL PARAMETERS ON SATELLITE ORBITS 




TABLE 1-4. SUMMARY' OF ORBIT ANALYSES FOR TIDAL PARAMETERS 


a 

e 

1 

k i 

k . s 

k n 

k n 6 

1 -30 

• 1 6 

33° 

• 312 

.0000 

-- 

-- 

1 .20 

.01 

49° 

• 250 

.0348 

-- 

m m 

1 .16 

.00 

90° 

• 342 

.0054 

•351 

.0074 

60.27 

• 055 

23-5° 

-- 

.0107 

-- 
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II. ORBIT ANALYSIS 


A. THE SELECTION OF ORBITS FOR THE DETERMINATION OF 
THE TESSERAL HARMONIC TERMS OF THE 
GRAVITATIONAL POTENTIAL 

When candidate orbits are selected for geodesy satellites, it is 
necessary to consider all possible ways in which measurable information 
can be obtained from such orbits. Existing satellites, of course, yield 
useful geodetic information. Drag-free orbits are therefore essential 
only in determining otherwise unknown and unobtainable data. Geodesy 
information obtainable from existing satellites is being considered, 
with results reported in Refs. 2-1 to 2-4. From these reports, it appears 
that the tesseral harmonic coefficients of degree 7-15 and order 3-10 
are the 72 pairs of coefficients which are least well known and probably 
of most current interest. Figure 2-1, taken from Ref. 2-4, indicates 
one group's opinion of the status of the tesseral harmonics. 

Improved values of the coefficients of order 11-16 have been 
measured by tracking the perturbations of existing satellites although 
the effect of atmospheric forces may cast some doubt on the validity of 
the data. 

Previous studies have discussed (and we have verified) the 
optimum orbits for measuring the unknown tesseral coefficients [Refs. 

2-7, 2-8], But it is of more practical importance to know the 
minimum number of drag-free orbits required to adequately determine all 
of these coefficients. 

In this Chapter, we review the pertinent perturbation theory and show 
how a drag-free orbit increases measurability by protecting long-period 
perturbations from distortions and allowing operation at lower perigee 
altitude. We conclude that 69 of the 72 pairs of coefficients of interest 
can be determined by Drag-Free Satellites in only four carefully selected 
orbits which are specified and discussed. 

1 . Analysis of Medium and Long-Period Perturbations . 

Medium-period perturbations will be defined as those with periods 

\ 

which are fractions of a day. Long-period perturbations are those with 
periods of several days. Perturbation analysis is begun by considering 
the gravitational potential expressed in the form 
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I 2 3 4 5 6 7 8 ? 10 II 12 13 14 15 

Order (m) 


FIG. 2t 1. STATUS OF KNOWLEDGE OF TESSERAL COEFFICIENTS 
ACCORDING TO REFERENCE 2r4 . 
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U 0 O i i 

V = r + e e r 

£ —2 m=0 p=0 q=-oo 


S R £mpq 


( 2 . 1 ) 


in which 


u / 3 e\^ 

R impq - a (a ) F imp (l) G £pq (e 


)f os | u 

(.sinj ( p- 


(i-m)even 


impq 

(i-m)odd 


( 2 . 2 ) 


and 


^mpq = ~ 2P)U + XI1 ~ 2P + q)M + m(a ~ 9 e “ (2 * 3) 


The quantities a, e, i, 0, oj, and M are the standard orbital elements, 

and u is the earth mass times the universal gravitation constant. 

J and are the amplitude and phase angle associated with the 

tesseral coefficients of degree i and order m. The angle Q is 

H 

the hour angle of Greenwich and a is the earth's equatorial radius. 

E 

The function F. (i), known as the inclination function, results from 
imp ’ 

the potential being rotated into the orbit plane. It is given by 
Kaula [Ref. 2-5] 



(2.4) 


where k is the integer part of (i~m)/2, t is summed from 0 to the 

lesser of p or k, and c is summed over all values for which the 

binomial coefficients are non-zero. The eccentricity function G. (e) 

ipq 

comes from converting orbit radius and true anomaly into a,e, and M. 
For (i - 2p + q) = 0, it is given by 


p«-l 


G. (e) = 

ipq 


i - i 


(l-e 2 /~2 ^0 \2d + l - 2p y i 



2p 


\ 2d+i-2p' (2.5) 


where p ' - 


p for p < i/2 


= (i - P) for p > i/2. 
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For (jj - 2p + q) £ 0, G. (e) is given by 

/ r ipq 

G . Pt ,< e > ” ( ‘ 1)|q| < 1 + p 2 /p l<,l 2 


k=0 


P 0 8 

ipqk ^pqk P 


2k 


where 


(2.6) 


P = 


fpqk 


(i + yrr^) ■ 

^(p'-i)^ ( _ 1) r fU - 2p' + 
2~t \ h-r / r.’ \ 28 


r=0 

h = k + q 7 , q ' > 0 
= k , q' ^ 0 ; 


and 


o = y (~ 2p ‘ ) — (u-z-2 

^ipqk . \ h - r/ r ! 00 1 ’ 


2 3 


h = J 


k , q' > 0 

k - a/ . a/ £ 0 . 


Here, 


p' = P, q / = q for p ^ &/2, 
p' = £~P, q' = -q for p > £/2. 


The instantaneous time rate of change of the orbital elements due 
to any one term in the gravitational potential is found from the 
Lagrange planetary equations 


a = 


JL 

an 3 m 


e = 


TlZ 


I2 1 


2 

na e 




3r 

3m 


3r\ 

3Tr) 


(2.7) 


COS 1 


w = - 


na2 JTe? s 


3 r Jl-e 2 ' 3 r 
~ + o — r- * 


in i 3i na^e 3e 
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1 = 


na 2 Jx-qV S in i 


. Sr Sr 

cos l-r— - -r— 

oco dn 


(2.7) 
Cont . 


n = 


l Sr 

na 2 Jx -e 2 * sin i Si 


M = n — 


i-e Sr 2 Sr 


na e Oe na- Sa 


where R = R. given in Eq. 2.2 and n is the mean orbital rate. 

i>mpq 

Approximate solutions to these equations are found by treating the 

it „ as linear functions of time [Ref. 2-5]. 

^mpq 

The medium-period perturbations of the in-track position of the 

satellite can be found for any given J . and from the closed-form 

£ m £m 

approximate solutions of Eq. 2.7 for the expression (A W + AM + A fi cos i). 
These perturbations, which fluctuate with a frequency of m times per 
day, are obtained by choosing particular combinations of the indices 
p and q such that (£ - 2p + q) = 0, i.e., the effect of the rapidly 

changing mean anomaly M is missing from the angle ^mpq' * or 
case, the amplitude of the in-track angular displacement, (AW + AM + 

Afi cos i) is 


(AW + AM + Afi cos i) 


max 


o 

iia J 
H E J flm 


na T (£-2p)w + m(fi-0 )] L 1 


7l-e 2 ' -(l- e^)_j F d£ + 2(U1)FG 


( 2 . 8 ) 


where the indices have been dropped from F. and G. . The rates 

£mp i!pq 

w and Cl are found by the usual first-order approximations 

• 3nJ 20Y\, . 2 0 

W = q - o [ 1 - 5 cos 1 J - 

4(l-e 2 ) a 


and 


ci = - 


3n J 2Q V 
2(l-e2) 2 a 2 


COS 1 . 
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Equation 2.8 is valid for any orbit, so this perturbation repre- 
sents a quantity which, if measurable, can be used to determine the 

coefficients J. and X . Additional perturbations to the cross-track 
im im 

and radial positions of an orbiting satellite also give a measure of 
the coefficients' magnitude, although the latter two are not as large 
as the in-track perturbation. In general, the medium-period fluctua- 
tions of the orbital elements are not large, and it is necessary to 
examine the long-period or resonant perturbations in order to obtain 
the desired coefficients more accurately. 

The closed-form linear solutions to the Lagrange planetary equa- 
tions all have the rate term 


i m Pq ~ ^ ~ 2p ^ W +(£ - 2p + q)M + m(ft - Q £ ) 


in the denominator. For instance, 

a E j 


2p a v £ 


£a = 


na 




FG(J - 2p + q) { 


£mpq 


cos 


sin 


(£ - m) even 

a 


111 , • 

) odd 


These solutions are valid except when this rate term is very small. If 
i]r , is small, the sine and cosine terms in the planetary equations 

v impq 

change very slowly, allowing large amplitude buildup. Then, the 

assumption of small perturbations leading to the closed-form solution is 

not valid. The condition if. ^ 0 occurs when (£ - 2p + q) m ~ m0 

Y tmpq . , E 

because co is small compared with M, and h is small compared 

with 9 ■ This is, of course, the condition of resonance for which the 
E 

orbital frequency is an integer multiple of the earth's rate. 

No general closed-form solution of the planetary equations, valid 
for the resonant case, is known at this time. Vagners [Ref. 2-6] 
obtained a solution for orbits of small eccentricity and ($-2p + q) = 1 
by combining the von Zeipel method with perturbation techniques of 
Hamiltonian mechanics. There are special cases where a pendulum solution 
for the osculating value of the nodal longitude of the mean satellite 
is quite accurate [Ref. 2-7]. In general, however, the exact solution 
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to the perturbations of the elements, due to those tesseral harmonic 
terms which are in resonance with the orbit, can be found only by 
numerical integration of Eqs. 2.7. The solutions depend on initial 
conditions and the degree of commensurability that the orbit has with 
the resonant condition. 

Consider the case where the number of revolutions per day equals 
the order m, that is, (£ - 2p + q) = 1. For this 


(q) = (1 - q)w + M + m 1 (fi - 0 E > 

= C L - qu) 


(2.9) 


where is a constant dependent upon the mean elements of the par- 

ticular orbit. Thus, for a given resonant orbit, there is a large set 
of beat frequencies ^ (q) corresponding to each q (assuming the 

inclination is not too close to the critical inclination. For the case 
of (£ - 2p + q) = 2 


i|t^ m (q) = (2 - q)W + 2M + m 2 (fi - 0 £ ) 

2 

= C g - qtb = 2C^ - qd) 

Hence, another large set of beat frequen&ies exists and, in fact, many 
other sets exist for ^ 0 or some multiple of to. Because of this 
spectrum of beat frequencies, it is theoretically possible to detect 
the effect of many tesseral harmonic coefficients from a given resonant 
orbit . 

The physical reason for the existence of the multiple beat fre- 
quency phenomenon is due to the combination of two effects. Orbital 
eccentricity causes a satellite to be more strongly attracted by the 
bulges due to tesseral coefficients at one part of the orbit (near 
perigee) than the other. In addition, the argument of perigee has the 
rate to which slowly moves the point of strongest attraction around 
the orbit. The tesseral harmonics have both zonal (i , latitude) and 
sectoral (m, longitude) dependence. The earth’s rotation makes a 
family of the tesseral harmonics (common m) resonate with the orbit. 
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The members of this family each have a different latitude dependence, 
however, and so the point of maximum attraction dwells for different 
durations and with different frequency near each. Thus, the long-period 
perturbations associated with a slightly-off resonant condition does not 
exactly coincide with the repeatability of the node passing through a 
wave length of the sectoral dependence; rather, the long-period effects 
are different, in general, for each latitude dependence. These effects 
show up in expanded powers of e and multiples of w. 


The rates ft, M, and u are not usually constant in the resonant 
situation, so, following Ref. 2-7> consider the acceleration of the 
elements. If the elements are designated x^ where i = 1, . . ., 6, 
then 


A. 


dt 


Sx. dx. ^x. 

"dt + ”3t 



All quantities but 
anomaly, 


M are of order J. . 

Am 


For acceleration of mean 
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jfcmpq 
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2 a jtmpq i m 
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( 2 . 10 ) 


^" 2P+ ^ F ^mp (i)G £ P q (e) 7m 


f-sin) , 1 ' ¥ °‘ 

{ > ^mpq 

l cosj ( £_ m) odd 


Although Eq. 2.10 is an approximation, it can be used to deter- 
mine the relative magnitudes of the perturbations to the in-track 
position of a satellite in a near-resonant orbit. Two remarks can be 
made about this equation: 

/ i 

(1) Because ~ (jig/a) , the effects of higher degree terms 

are attenuated. 
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2. G £pa^ e ^ is P r °P° rtional to e . Thus, higher values of I Q ! 
tena to decrease Mf,mpq* This is especially true for small 
eccentricity. The highest possible eccentricity is therefore 
desirable. 

For a resonant orbit, the number of pairs (J. , X . ) of unknowns which 

& m I m 

can be determined is equal to the number of distinct beat frequencies 
which can be measured from the in-track perturbation. This is true, 
provided that the satellite remains in. the orbit for at least one-half 
the period of the smallest required beat frequency. Thus, it is important 
that the orbit be "tuned" closely enough to exact commensurability so 
that the resonant effect is a large one, and yet preferably not so close 
that one beat period is too long to observe in the satellite lifetime. 

To determine which resonant orbits should be used to obtain particular 

(J. , X. ), the magnitude of M„ is computed (using hypothetical 

f.m 1 m rmpq 

values for J. ) and examined for measurability. The eccentricity is 
xj m 

assigned a value which seems very comfortable for a Drag-Free Satellite 
and the dependence upon inclination is investigated by letting inclina- 
tion range from 0° to 90°. 


2. The Se lection of Resonant Orbits 

The above analysis suggests that the following procedures should 
be used in selecting resonant orbits for determination of poorly known 
tesseral harmonic coefficients. 

(a) Determine which tesseral harmonic terms can be adequately 
measured by observing the medium-period perturbations of 
existing satellites. 

(b) For important terms not obtainable from existing satellites, 
investigate the possibility of determination from satellites 
placed in resonant orbits. The order (m) of the terms deter- 
mines the semi-major axis of the orbit required. The inclina- 
tion of the orbit should be chosen so the resulting accelera- 
tions in mean anomaly (M) for each term of order m is large 
enough to produce a detectable effect. The total effect can 
be determined by computing M vs i for the various combina- 
tions of (i , m, p, q) of interest. For example, suppose 

it is determined from a series of M vs i plots that an 

inclination exists such that in-track perturbations are 

measurable for q in the range -3 to +3. Then seven 

pairs (J, , X, ) can be determined from this orbit. If 

£m’ j|m 
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\ 


there are ten pairs of coefficients required with this 
particular order m, then at least two orbits will be 
required for this particular semi -major axis. 

(c) This procedure can be repeated for each m in the range of 

orders of the terms sought. If a range of 

acceptable inclinations exists for any two adjacent values of 
m, then satellites can be put into both these orbits, using 
the same booster, without a plane-change maneuver. 

(d) If a desired pair of coefficients (JgmAiJni) can be only 
marginally determined from a particular resonant orbit, a 
comparison must be made of the relative merits of establishing 
that orbit or trying to determine the coefficients from the 
medium-period perturbations on another orbit. 

As indicated before,, coefficients of degree 7-15 and order 3-10 
are of primary interest. It is feasible to determine all these coeffi- 
cients from four satellites which have orbital frequencies of 3, 4, 5, 
and 7 times per day. 


3. Computation of The Amplitude of Mg m p q 

I m| vs i computations have been made for different coefficients 
1 1 max 

in both Refs. 2-7 and 2-9. In these studies, however, the value of q 
was limited to 0 and + 1. Here, because we are looking for a measure 
of the perturbations due to many beat frequencies, these computations 
need to be expanded. 


The hypothetical value of 

J. used 

£m 

in 

the 

computation of the 

amplitude of is 

iJmpq 



1/2 


J ^m 

a - 

m).’ (4£ + 

2) 
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(£ + m): 
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£m 
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where 

• £™ 

approximat ion 


is 


the normalized coefficient. 


Kaula [Ref. 2-5] gives the 



and Allan [Ref. 2-9] uses 


J2 X 10 


-5 


f 


( 2 . 12 ) 
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10 


-10 


£m 


12.2 X 10 
m + l) 2 (2i + 3) 


(0.93) 


2£+3 


1/2 


which are nearly equivalent. Although published estimates of actual 

(J. , A. ) through (15, 15) exist these only tend to confuse this 
15m jjm 

analysis, so Eq. 2.12 is used instead. 

To determine coefficients of degree 7-15 from a single orbit, at 
least nine beat frequencies must be detectable. To get an idea of how 
the perturbation is attenuated by increasing |q| , a series of computer 
runs of |M| vs i was made for m = 7, £ ranging from 7 through 15, 

and q ranging from -4 to +6. The altitude of perigee was chosen as 
450 km which is easily within the capability of a long-life Drag-Free 
Satellite, The semi-major axis was established as that required for a 
satellite to orbit the earth seven times per day. The output of these 
runs consists of a series of plots of | M j (denoted by |mDD0t|) vs i 
and is found in Appendix A. Here, the index s is the number of orbits 
per day. The quantity (£ - 2p + q) in Eq. 2.10 is replaced by (m/s). 

Referring to Appendix A, it is seen that, in general, the magnitude 
of j M | does indeed decrease with increased £ and | q| . It was desired 
to determine how large the in-track resonant orbit perturbations are 
due to the relevant coefficients of order m = s. It also was necessary 
to determine if measurable perturbations exist corresponding to each 
beat frequency. Consequently, the index q - 0 was assigned to the 
(|,p) = (15,7) combination and other q's were assigned to other com- 
binations of lower degree £. 

Another series of computations of |m | was made for resonant orbits 
with s = m ranging from 3-6 and 8-12 with altitude of perigee at 
450 km. The combinations of ($,P,q) investigated in particular were 
(15, 7, 0), (14, 7, 1), (13, 7, 2), (12, 5, -1), (11, 4, -2), (10, 6, 3), 
(9, 6, 4), (8, 2, -3), and (7, 1, -4) although some other combinations 
were included. The plots from these runs are presented in Appendix B. 

By examining the plots of Appendices A and B, it can be seen that 
the general level of |m| over the range of inclinations increases with 
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increasing order m. This is primarily due to the decrease in semi- 
major axis as m is increased with perigee altitude held constant. It 
can also be seen that there are several places where |m| sharply decreases 
and then increases as inclination is changed. These places correspond 
to points where the inclination function F(i) changes signs. Such 
crossover inclinations should be avoided. If one removes from consider- 
ation, the near neighborhood of critical inclination, those inclinations 
corresponding to crossover points, and those with M less than what can 
be regarded as a minimum detectable level, then what remains are those 
ranges of i constituting possible inclinations for a geodesy satel- 
lite. 


If the minimum beat period is 60 days, an acceleration with maximum 
-5 . 2 

value of 10 deg/day will result roughly in an intrack perturba- 
tion of amplitude between 100 and 1,000 meters. Thus, one can consider 

an |M | of 10 deg/day or larger as constituting a substantially 

| -6-5,2 

measurable effect, |M | between 10 and 10 deg/day as being margin- 
ally detectable, and |m| less than 10 ^ deg/day^ as being undetectable 
Applying this criterion to the data of Appendix B, one can see that the 
J , J , J , J , and J terms for the combinations of 

It) j u It j *J Idj O i j U IO , T 

(p,q) chosen are only marginally detectable. (In fact, the general 

magnitude of the acceleration of J. _ „ is so small that unless close 

15 , o 

tuning can be achieved, it can be as accurately predicted from medium 
period effects on low satellites.) The rest of the coefficient combi- 
nations tend to have large ranges of quite detectable accelerations. 


As an example of the mean anomaly accelerations |m | , for the 
distinct combinations of (Jh, m, p, q) investigated, Table 2-1 was prepared 
for various values of m. Two values of inclination are shown for the 
orbits with s = 3-6. For the inclinations shown, only J is not 

lo ^ o 

detectable, J UJ , J 13>3> J 12>3 . J v , 3 . J is, 4 ’ J 7,4> ond J 7,6 are 
marginally detectable, and the rest are easily detectable by the cri- 
terion given above. This is dependent, of course, upon the particular 
combinations of (p,q) chosen. 


Consider now the possibility of determining coefficients from orbits 
where m = 2s, 3s, and 4s, i.e., conditions of greater than one-day 
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TABLE 2-1 

MAXIMUM ACCELERATION OF MEAN ANOMALY, | !l| X 10 5 DEG 
FOR GIVEN INCLINATION AND INDICES L, M, P, Q WITH S 


/ 



10 3 









commensurability. Since 


M 


i m Pq 


3(i ,a. 

S 

a 


,a . & , m v 

(i) (:) 


F. (i) 

£mp 


G„ (e) 
iPQ 


£ra’ 


then for fixed (JJ, m, p) , the attenuation of |m| for the higher orbit 
will be 



(2.13) 


Here, the indices 1 and 2 represent those quantities corresponding to 
the lower and higher orbit respectively. The value of q corres- 
ponding to G (e) is equal to (q + S /s ) where q. corresponds to 

/ i 1 4 i 

(£ + 2p + q^) = m/s^. Tables 2-2, 2-3, and 2-4 show ratios of | | / | | 

for m = 6, 9, and 12 in the s = 3 orbit, m = 8 and 12 in the s = 4 

orbit, and m = 10 in the s = 5 orbit. Table 2-5 shows the resulting 

values of |?»ijfor inclinations of 58° for s = 3, 50° for s = 4, and 

52° for s = 5 orbits. It can be seen that except for J . all 

15 , 6 

overtone coefficients are at least marginally detectable in the s = 3 

orbit. For s = 4, only the J term is marginally detectable with 

12,0 

the rest being quite detectable by the above criteria. For s = 5, 
all m = 10 terms are quite detectable. It must be remarked that no 
special effort was made to maximize the effect of the overtone terms 
in selecting the inclinations used for Table 2-5. With further effort, 
there can undoubtedly be improvement. The point is that a great deal 
of information can be obtained from the s = 3, 4, and 5 resonant 
orbits in addition to what is obtainable from the fundamental beat 
frequencies. In fact, all but three of the 72 pairs of the desired 
coefficients can be found from those orbits. It also must be emphasized 
that the choice of cutoff points between what are referred to as undetect- 
able, marginal, and quite detectable perturbations is strongly dependent 
upon uncertainties due to tracking and perturbing surface forces on the 
satellite. For a Drag-Free Satellite, the effect of surface forces 
can be attenuated to the point where tracking uncertainty is the major 
limiting factor. 
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TABLE 2-2. ATTENUATION OF THE m = 6, 9, and 12 COEFFICIENT ACCELERATIONS IN THE S = 3 ORBIT. 
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TABLE 2-3. ATTENUATION OF THE m = 8 and 12 COEFFICIENT ACCELERATIONS IN THE s = 4 ORBIT. 
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COMPARISON OF MAXIMUM ACCELERATIONS OF THE MEAN ANOMALY 
5 2 

|m| X 10 DEG/DAY , DUE TO HIGHER OVERTONE TERMS 
(s = m), FOR s = 3, 4, and 5 REVS/DAY. 



fH 




V 







o 













rH 

00 

' CO 

oo 


00 


Cl 




II 

O' 

CM 

o 

CO 

co 

co 


CO 




e 

ii 

• 

• 

• 

• 

• 


• 





CM 

CO 

O 

h- 

co 

o 


CO 





Of 


rH 

pH 

CM 

CO 


rH 





o 

CM 

rH 

o 


CO 

■O’ 


o 

co 

CO 


CM 

CM 


o 


TT 

CO 


CD 

in 

rH 

to 

m 

• 

• 

• 

• 

• 

• 


• 

• 

• 

II 

H 

rH 

CM 

o 

ci 

o 



rH 


m 

G 




rH 


rH 

CM 







o 

rH 

CM 

rH 

CM 

CO 

CO 


CO 



o 




1 * 

1 


1 

1 

1 


in 

it 

w 

a, 



h* 

in 


co 

CO 

CM 

CM 

m 


in 


CO 

CM 

rH 

o 

o 

Cl 

00 

N 



rH 

rH 

rH 

rH 

rH 

iH 

rH 





m 


CO 

co 

CM 

CO 

co 

in 

CM 

b- 

CO 

rH 

b- 

CM 

• 

• 

• 

• 

• 


• 1 1 



CO 


rH 

CM 

CM 


t" CD 

cn <£> 


ci o 
co in 


oo co 

co m o 

II II 
G -rl 


m co cm 
in in 


oo co 

cm o 


CM 

rH 

CM 

CO 

CO 


CO 


1 

1 


1 

1 

1 

b- 

1/5 

TT 

co 

CO 

CM 

CM 

00 

CM 

rH 

o 

o 

05 

00 

rH 

rH 

rH 

rH 

rH 
















4 . Additional Advantage of the Drag-Free Satellite for Orbits of 

Small Major Axes . — — 

To determine the coefficients of the harmonic expansion of the 
earth's gravitational potential from satellite perturbations, it is 
important that these perturbations not be obscured by other forces such 
as those due to radiation pressure and atmospheric pressure. Because 
these forces are always present to some degree, their effective removal 
by use of a Drag-Free Geodetic Satellite will always produce more 
confidence in the validity of the tracking data. The situation in 
which greatest improvement would be obtained is with orbits of low 
perigee altitude. 

In the literature [Refs. 2-1 through 2-4], it is reported that 
there are a number of resonant satellites in the range s = 11-15 and 
that the coefficients with these orders can be reasonably well deter- 
mined. If one considers a circular orbit with a frequency of 14 rev/day, 

q 

because of the e 1 factor in the G(e) function, only the beat fre- 
quency corresponding to q = 0 is detectable. One cannot assume that 
for £ > 14, the tesseral coefficients are so small that their effect 
is negligible if Eq. 2.11 and 2.12 are valid. In other words, for a 

circular orbit, the J J and J, „ , . coefficients will all 

15,14 17,14 19,14 

significantly affect the q = 0 perturbation. Thus, one must make the 
orbit as eccentric as possible to obtain these coefficients by increas- 
ing the number of detectable beat frequencies. A Drag-Free Satellite 
can contribute significantly here by allowing the low perigee altitudes 
necessary for increasing eccentricity to an acceptable value. 


To illustrate the effect of an increase in eccentricity, Table 2-6 
was prepared showing the change in the magnitude of the ^mpq^ 6 ^ 
function for the s = 14 orbit. The values of degree £ ranged from 
14-20 and q from -3 to +3 which would produce 7 beat frequencies. 
The eccentricities used correspond to perigee altitudes of 350, 500, 

650, and 800 km. Appendix C contains the [m | vs i plots corresponding 
to this range of indices for an altitude of 350 km. To determine the 
| M | for any other perigee altitude, multiply by G(e)/G(0. 0748) . 

For example, when i = 50, the following results are obtained in a 
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TABLE 2-6. COMPILATION OF THE G(e) FUNCTION FOB THE s = 14 ORBIT CORRESPONDING TO PERIGEE 
ALTITUDES OF 350, 500, 650, and 800 km. 



comparison of the 350 and 800 km perigee altitude cases. 





TABLE 2-6 


COMPARISON OF MEAN ANOMALY ACCELERATIONS FOR 


FOR 

PERIGE 

ES OF 350 and 

800 KM 




deg/day2 

i 

P 

q 

|^35o! X 10 5 

|^80ol * 10 

14 

5 

-3 

2.752 

0.013 

15 

6 

-2 

32.754 

0.850 

16 

9 

3 

1.240 

0.006 

17 

9 

2 

24. 748 

0.622 

18 

8 

-1 

196.01 

26.38 

i 19 

9 

0 

332.51 

209.20 

1 ' 20 
1 

10 

1 

171.72 

■ 

21.82 

. 


It is evident that the eccentricity effect is significant for j q | > 2 
The perturbations due to coefficients of degree £ > 15 seem to be 
substantial . 
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III. SATELLITE DESIGN 


Since 1963, Stanford University has been developing drag-free tech- 
nology for application to several possible missions (e.g., an aeronomy 
mission, and an unsupported gyroscope mission). Throughout this work, 
we have become convinced that it is possible to build the basic mechanism 
(proof mass, pickoff, processing electronics, and thrusting system) 
with standard state-of-the-art flight hardware and techniques. The 
interface of the drag-free device with the rest of the satellite can 
be simple or complicated depending on the requirements for the other 
aspects of the mission (other experiments, physical configuration con- 
straints, kind of telemetry desired from the drag-free mechanism, etc.). 
Reference 3-1, for example, proposes a very simple complete Drag-Free 
Satellite for use in a low altitude aeronomy mission. 

The discussion in this part of the report has therefore been limited 
to the feasibility of applying the drag-free principle to goedesy missions. 
In particular, it is concerned with any special requirements arising 
through application to geodesy. 

The drag-free device very effectively cancels the disturbing accel- 
erations due to surface forces such as atmospheric drag and solar radia- 
tion pressure. The largest remaining non-geodetic disturbing acceleration 
is that due to the mass attraction of the satellite itself on the proof 

mass. This acceleration can be as low as 10 ^g or as high as perhaps 
“8 

10 g depending on the care with which one manages the mass distribution 
within the satellite (particularly the masses within say, 10 cm of the 
proof mass). 

The principal component of this mass attraction is due to the fact 
that the satellite mass center does not necessarily coincide with a point 
of zero mass attraction. Since the proof mass is nominally at the 
satellite mass center and the points of zero mass attraction are fixed 
points in the satellite, spinning the satellite about an axis normal to 
the orbit plane will tend to average the orbit plane components of this 
force. This will reduce the intrack disturbing acceleration by about 
2 orders of magnitude below that attainable without spin. 


- 50 - 



Sections A and B which follow , present solutions to problems 
which arise for a spinning Drag-Free Geodetic Satellite. First, active 
attitude control to maintain the spin axis normal to the orbit plane; 
then, phenomena in the translational control uniquely associated with 
spin and the system mechanization are discussed. 

It is not necessary, of course, that a Drag-Free Geodesy Satellite 
be spin stabilized. Quite useful geodetic information could be obtained 
with a gravity stabilized Drag-Free Satellite, such as GEOS-C with a 
modular, add-on drag-free package. With such a satellite, there is the 
question of whether or not the translational control system could couple 
into attitude motion (through misalignment errors in the pickoff and 
thrustors) in such a way as to cause attitude instability. This topic 
is discussed in the final section of this Chapter. 

A. MAGNETIC ATTITUDE CONTROL OF A SPINNING DRAG-FREE GEODESY SATELLITE 

As discussed in Chapters I and II, intract perturbation of a 
Drag-Free Satellite is the effect which contains the most significant 
geodetic information. It is desirable to minimize as much as possible, 
the effect of all other intrack disturbances from these measurements. 

This can be done most easily with a spinning Drag-Free Satellite with 
spin axis normal to the orbit plane. The spinning motion allows very 
effective averaging in the orbit plane of the major components of mass 
attraction of the satellite mass on the proof mass. 

Because many disturbance torques are present which can drive the 
spin axis from the orbit plane normal and change the satellite's spin 
speed, an attitude control system must be provided to correct for these. 
An excellent method of providing this attitude control for the orbits 
of interest is to make use of the earth's magnetic field. By creating 
a magnetic dipole moment in the satellite, a torque is produced on the 
satellite as this moment tries to align with the earth's magnetic field. 
By controlling the direction of the dipole moment, one has the ability 
to control both the pointing direction and spin speed of the satellite. 

The basic idea of magnetic attitude control is not new and has been 
studied extensively and used in several satellite systems. In the 
following sections, a new method of three-axis magnetic attitude control 
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is developed which can rely on a single magnetic coil. First, the 
linearized equations of attitude motion are developed and the effect of 
a nutation damper is added. Then, a brief review of disturbance torques 
acting on the satellite is presented. Following this is a discussion 
of the use of a state estimator driven by horizon sensors to determine 
the attitude errors. With this information, the magnetic attitude control 
system is developed for a satellite with mass symmetry about the spin 
axis. This system has three modes so that both pointing control and 
spin control can be achieved. A procedure for determining control gains 
which will provide an average optimum performance is presented. Stability 
is demonstrated by using Lyapunov and averaging techniques. Then, a 
brief summary is made of the analog and digital simulations of this 
attitude control system. 

1. The Attitude Dynamics of Rigid Spinning Spacecraft 

In this part, equations representing the dynamics of the spinning 
spacecraft with respect to relevant coordinate systems are developed as 
a means of notation and coordinate system definition. 

a. Spacecraft Kinetics . 

If it is assumed that the control axes (the axes about which 
control torques can be applied) are the principal inertia axes for the 
spacecraft mass center, then the kinetic equations describing the 
vehicle's attitude motion are the well-known Euler equations 
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Here, (I„„ , I,,,,, I„„) are the principal moments of inertia. The 

' ax yy 

vector ujB"! is the instantaneous angular velocity of the body with 

respect to an inertial reference frame. The quantities (to , to , to ) 

x y z 

are the measure numbers in the body-fixed frame aligned with the principal 
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axes 


. The quantities (T , T , T ) are the measure numbers in the same 

x y z 

body-fixed frame of the external torques acting on the satellite. 

b . Coordinate Systems and Spacecraft Kinemat ic s . 

The coordinate systems of interest are shown in Figs. 3-1 
to 3-3. The primary inertial reference frame is defined by the x^. 
axis pointing to the Vernal equinox, the z^ axis pointing along the 
earth's spin axis, and y^ completing the right-hand orthogonal set. 

Figure 3-2 shows the "local" reference frame with x along the 

L 

direction of the vector R from the geocenter to the satellite. The 

axis z is normal to the orbit plane and along the orbit's angular 
L 

momentum vector, and y completes the right-hand orthogonal set. 

The angles i, ft, g, and f have the usual definitions of inclination, 
right ascension of ascending node, argument of perigee, and true anomaly. 
The local frame will be used as the reference set (R) later, in the 
study of controlling alignment of the spin axis with the normal to the 
orbit plane. 

c . Linearized Equations of Motion 

When it is desired to keep the spin axis aligned normal to the 
orbit plane, the natural axes about which to apply control torques are 
the local axes (L) of Figure 3-2. These correspond to the reference 
axes (R) of Figure 3-1. However, the closest possible positions of the 
principal axes of the satellites are the intermediate axes labeled (P) 
in Figure 3-1. These axes correspond to the orientation of the body- 
fixed axes at the time when the angle equals zero. 

Assume that the angles 0 and Q (corresponding to yaw and roll 
angles about x R and y p respectively in Fig, 3-1) are small enough so 
that the approximations 

sin 0 = 0 , 

sin 0 = 0 , 

cos O = cos 0=1 , 

can be made. 
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FIG. 3-2. ORIENTATION OF THE LOCAL (L) 
REFERENCE FRAME 



FIG. 3-3. GEOMETRY OF THE ANGLES AND 
ANGULAR RATES OF THE WOBBLE 
PROCESS. 



• » ♦ 

Assume also that w = * = constant » g + f. Then, for a satellite 

z 

with mass symmetry about the spin axis (i.e., I = I ) . Ea 3.1 

xx yy 

can be written in the constant matrix form 
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(3.2) 


In this equation, a and a are the x and y measure numbers of 

x y P ‘ P 

the angular velocity of the body with respect to the inertial frame, 
and T x and T v are the x and y of the torques applied to the 
body (normalized by I ). D is the ratio of the moments of inertia 

XX 

times the spin rate, 


D = I co /i 

zz z' xx 


The constant n corresponding to mean orbital rate, has also been 
substituted for g + f, which is time varying for an elliptic orbit. 

d . Addition of a Wobble Damper . 

Figure 3-3 depicts the geometry of the angles and angular rates 

of the wobble process. The body axes x r y , z are rotated from an 

B B B 

inertially fixed set X, Y, Z through the classical Euler angles 

tj , and Q. The Z axis is aligned with the total angular momentum 
«. • • • 

vector H. The rates § and T) are the inertial precession and nuta- 
tion. The nutation angle V is described by 

cos tj = I w /H . (3.3) 

1 zz z' 

The term "wobble damping" is defined as the driving of tj to zero. 

If the nutation angle tj is small, and u is held approximately 
constant, then the so-called "energy sink approximation” technique 
yields the relationship tj = -drj where d is the time-constant of 
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the wobble damping process. 


After a bit of algebra, it can be shown that the effect of 
the damper is the addition of two diagonal terms in the 4X4 system 
matrix of Eq. 3.2. Equation 3.2 thus becomes, 



— — 


— — 

-d -D 0 0 


a 

X 

D -d 0 0 


a y 

1 0 0 n 


♦ 

0 1 -n 0 


e 



(3.4) 


Equation 3.4 will be used throughout (this) section, A, to determine the 
desired control system for a spinning symmetric geodesy satellite. 

The exact equations of a mechanical wobble damper are generally more 
complex than has been assumed by the energy-sink procedure used here 
and require evaluation of a higher-order system [Refs. 3-2, 3-3, 
and 3-4] . However, most references on the subject indicate that the 
energy-sink approximation gives results reasonably close to actual 
behavior, so the method will be retained rather than increasing the 
order of the system equations which would be necessary for a closer 
investigation of the effect of a particular type of damper. 


2. Disturbance Torques 

Disturbance torques acting on the satellite can be broken into 
three categories. Two, referred to as "inertially-f ixed" and "body-fixed" 
torques, are those which tend to move the spin axis from the reference 
and are about the vehicle's lateral axes. The third type of torque 
is about the spin axis and tends to change the rate of spin speed. 

Short discussions of the following important inertially-f ixed 
torques will be presented here 

a. atmospheric 

b. radiation pressure 

c. misalignment of translation control jets 

d. magnetic effects 

e. reference-frame kinematics. 
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Full developments can be found in the indicated references 


a. Atmospheric torques 

Consider the atmospheric torque acting on a cylindrical satel- 
lite with uniform wall material. If 6z is the distance from the 
mass center alogg the axis of symmetry to the geometric center, then, 
from the analysis of Ref. 3-5, it can be shown that the inertially 
fixed torque magnitude acting on the satellite is 


T 

aero 


(T end + 


P cy! )6z 


(3.5) 


where, for diffuse reflection 

P cyl = ^ r P V i 2 sin t, l [sin t, l + 1T/6 (V^)] 

2 o 

T = TTr pV. sin q. cos q. 

end H i *i '*i 

and 

l = cylinder length 

r = cylinder radius 

p = atmospheric density 

= relative speed of incoming molecules 

V = relative speed of outgoing molecules 
r 

tj = angle between incoming velocity and the wall surface. 

The ratio V /v — Jl - Cl ' , where a is the accommodation coefficient 
r' i 

for the particular surface. References 3-5 and 3-6 present accommodation 
coefficients ranging from 0.3 to 0,95 depending on the wall material 
and gaseous medium involved. The relative velocity is found by 


where V 
w 

e 

R 

K 


V. = -V + KU) X R 
l e 

inertial velocity of the spacecraft 
earth spin velocity 
radius vector from the geocenter 
wind constant. 


(3.6) 
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I 


The angular velocity of the upper atmosphere has been determined by 
examining changes in inclinations of various satellites by King-Hele 
[Ref, 3-7]. He found K = 1.46 for nine satellites at heights of 
200 to 300 km. The atmospheric model used here is developed in 
Refs. 3-8 and 3-9. 

b. Radiation pressure torques . 

Direct solar radiation pressure in the vicinity of the earth is 

Ra^ = 4.66 X 10 5 dyne/cm 2 . (3.7) 

Earth emitted radiation pressure is 

Ra == 7.53 X 10~ 6 dyne/cm 2 . (3.8) 

. *3 

By numerical integration of the solar energy reflected from the earth 
to a satellite, the following emperical relationships have been found 
for the pressure due to this reflected energy. 

—5 -4 2 

Ra = 1.80 X 10 exp(-3 X 10 h) cos {3 dyne/cm , (3,9a) 

O 

in which h is the satellite altitude in km, and p the angle between 
the earth-satellite radius vector and the earth-sun line. Ra is 

O 

assumed zero for (3 > T! /2. The reflected radiation vector is at an 
angle of (p + v) rad from the earth-sun line, where 

V i fiO/ff) 2 * 4 rad (3.9b) 

and 

f 4 4.89 - h(5, 82 X 10~ 4 ) rad. (3.9c) 

According to Evans [Ref, 3-1 0] , the pressure and shear stress 
components due to a radiation source vector Ra striking a wall at an 
angle Q range from 
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p = (Ra) sin £ (sin £ + 2/3 p) 

t = (Ra) sin £ cos £ 
for diffuse reflection to 

2 

P = (Ra)(l + p) sin £ 

T = (Ra) (1 - p) sin £' cos £ 


(3.10) 


(3.11) 


for spectral reflection. Here, p is the surface reflectivity. Because 
the radiation pressures can be formulated as vectors, the radiation 
disturbance torque evaluation is done in exactly the same fashion as 
for aerodynamic torques. 

c. Torques due to translational control thrustors . 

There are two apparent ways in which translation control jets 
could cause disturbance torques to the satellite attitude. If the line of 
action of a translational control thrustor does not pass through the satel- 
lite center of mass, the system will produce torques which could fleet both 
the pointing accuracy and the spin speed. In addition, a leak in a pneu- 
matic system either at a joint or in a valve could also produce torques. 

d . Magnetic torques . 

Magnetic disturbance torques are primarily caused by current 
loops in the spacecraft and materials subject to permanent or induced 
magnetism. The instantaneous torque is the vector cross product of 
the spacecraft's effective dipole moment M^ and the magnetic induction 
of the local field B or 


T = M X B . (3.12) 

m s 

Here, M is considered to include all but the dipole generated by the 
s 

control coils. For a spinning satellite, the inertially fixed torque 
will equal the product of the dipole component along the spin axis 
and the lateral component of the magnetic field. 
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Because the satellite spins with respect to the magnetic field 
vector, torques due to the induced currents (eddy currents) and the 
irreversible magnetism of permeable materials (hysteresis effects) 
must also be considered. Smith [Ref. 3-11] has developed a theoretical 
expression for determining the eddy current torques in rotating shells 
which is used in this analysis. 

The magnetic field assumed for the analysis of these disturbance 
torques is modeled in Ref. 3-12. Coefficients used in this analytical 
model are presented in Ref. 3-13. 


i , the 

(3.13) 


e. Reference frame kinematics. 

For a satellite in an earth orbit with inclination 
instantaneous orbit precession rate is, to first order, 


n = - 


2 2 
3^i J R cos i sin o 

Ct © 

HR 3 


Here, p is the universal gravitation constant times the earth mass, a, 
is the sum of the orbit’s argument of perigee and the true anomaly, H 
is the orbit's angular momentum with respect to the earth, and J 
is the first harmonic term in the expansion of the earth's potential. 

If the spin axis of the satellite is to be maintained normal to the 
orbit plane, a torque, "kipematic disturbance torque," must be exerted 
on the satellite to precess the spin axis as the orbit normal precesses. 
A simple consideration of required rate of change of satellite angular 
momentum yields the torque expression 


T 


K 


0, sin i 


I co 
zz z 


K 


(3.14) 


where y is a unit vector along the line of nodes. 

K 

f . Total torques 

To determine the total disturbance torque acting on any 
particular cylindrical spinning satellite, a digital computer program 
was written to evaluate these five torques. This program computes the 
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torques as a function of the parameters which change with orbital posi- 
tion of the satellite in its nominal orientation. Examples of iner- 
tially fixed yaw and roll torques were computed using this program and 
are presented as functions of orbital position in Figs. 3-4 and 3-5. 

For these examples, the orbit is inclined at 45° and the orbital fre- 
quency is 15 times/day. Perigee occurs at 300 km over the equator 
which is positioned in the center of the atmospheric bulge. The dimen- 
sions and other parameters used for the cylindrical spacecraft of this 
example are 

length i = 0.5 m, 

radius r = 0.5 m and 0.75 mt (two cases), 

skin thickness = 0.1 cm of aluminum, 

spin speed = 1 rad/sec, 

accommodation coefficients = 0.64 and 0.84 (two cases), 

decimeter flux index = 300 and 100 (two cases), 

reflectivity = 0.0 and 1.0 (two cases), 

2 

spin axis dipole = 0.6 amp m , 

c.m. offset from c.p. = 0.5 cm, 

jet misalignment moment arm = -0.5 cm. 

Such plots are useful for illustrating the relative magnitudes 
and characteristics of the various disturbance torques. 

3. State Estimation by Filtering Horizon Sensor Data 

By proper placement of a pair of infrared horizon sensors (bolometers), 
it is possible to measure directly the roll error, 0. Because of the 
structure of the satellite's state equations, the system is "observable." 
That is, a filter can be constructed which will produce estimates of 
the other state variables from the roll error measurements. To account 
for driving noise (disturbance torques) and measurement noise, an 
optimum steady-state (Kalman) filter is developed which gives "best" 
estimates of these states. This system is developed as follows. 

a . Horizon-sensor determination of roll angle, spin speed, and 
orbital rate. 

Roll- error measurement can be provided using the output of 
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two infrared bolometers arranged with their optical axes in a "Vee" 
configuration as shown in Fig. 3-6. The optical axes lie in a plane 
containing the vehicle spin axis. As the satellite spins, the optical 
axes sweep out two conical surfaces in space. The sensors produce 
signals related to the change in received radiant energy as optical 
axes sweep from cold outer space through the warmer infrared earth, 
and back to space again after each spin revolution. 

The intersection of these sensor paths with the earth and the 
corresponding sensor output is shown in Fig. 3-7. When the spin axis 
is normal to the earth-satellite radius vector, the sensor pulse out- 
puts will have the same width and occur at the same time for an ideal 
spherical earth. If a roll error exists, the relationship among pulses 
will be as shown. Errors in yaw cannot be instantaneously detected. 
However, since the satellite spin axis is approximately fixed in inertial 
space, a pointing error in yaw at any particular instant will become 
an error in roll 90° later in the orbit due to the rotation of the 
local reference frame. 

The geometry of the horizon sensor scheme is dipicted in Fig. 


3-8 where 



e 

= 

roll error 

6 

= 

half-vee angle 

a 

= 

90° - 6 

P 

. = 

sin ^ (Re/R) 

t i 

= 

time between horizon pulses for right 

• 


, , 



local spin rate = co - f - g 

z 

y 

= 

t^/2 . 


The earth sweep time for the right sensor, t^, can be written 


t 


1 


2 -1 
t- cos 

* 


cos 5 + sin 9 sin 6 
cos Q cos 6 


(3.15) 


Equation 3.15 was used to compute the output characteristics 
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FIG. 3-6. HORIZON SENSOR "VEE" CONFIGURATION 



FIG. 3-7, HORIZON SENSOR OUTPUT WITH AND 
WITHOUT ROLL ERRORS. 
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FIG. 3-8, GEOMETRY OF THE HORIZON 
SENSOR SCHEME. 
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of sensors with half-vee angles of 3° , 6° , and 9° for an elliptic orbit 

having 3 revs/day. Define At and T as 
' avg 



where t is the time between horizon pulses for the left sensor. 

2 . 

and T = (t + t )/2. Then Fig. 3-9 shows roll-angle error as a 
avg 2 r 

function of . At for a 9° half-vee angle for various values of 
satellite radius. The roll error is found by multiplying At by a 
gain (slope of curves such as those in Fig. 3-9) which is radius de- 
pendent. This gain, as a function of T , is shown in Fig. 3-10 

avg’ 

for the three sensor angles investigated here. The quantity T aV g 
average remains fairly constant at a given altitude for small values 
of roll error. As can be seen from Eq. 3.15, the sensors will produce 
a larger error signal with a larger half-vee angle. However, a 
trade-off must be made because larger sensor angles will cause the 
sensor to miss the earth for high altitudes of highly eccentric orbits. 


B.ecause the horizon sensors produce a pulse train, with T a 
a function of satellite altitude, this output can also be used to give 
a direct measurement of the vehicle's spin speed and the current orbital 
rate, f + g, of the vehicle. 

b . Kalman filter for state estimation . 

The linearized state equations of the symmetric spinning satel- 
lite was summarized in matrix form in Eq. 3.4. This equation, as it 
stands, is observable from a measurement of the roll error Q. That 
is, with a signal proportional to the roll error Q, a state observer 
can be constructed to estimate values of the yaw error <!> and the two 

inertial rates a and Ct • 
x y 

It would be nice if the disturbance torques could also be 
treated as states and estimated as part of the observer's function. 
However, a constant torque about the x-axis (yaw axis) (e.g., aerodynamic 
torque) is not observable from roll error measurements of a symmetric 
vehicle. Therefore, for the time being, it will be assumed that the 
disturbance torques can be treated as white driving noise. If it is 
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assumed that the driving noise and measurement noise are stationary, 
v/hite, uncorrelated processes, a Kalman filter [Ref. 3-14] can be mechan- 
ized to determine the linear least squares estimates of a > Ct > anc * 

x y 

4 >. 

The state and observation equations are expressed in general 

form as 


x = Fx + Gu + Lv 
y = Hx + w 


(3.16) 


where u is the normalized control input, and v and w are white 
noise with covariances 


E(v(t)v T (r) ) = Q6 (t - t) 

E{w(t)w T (T>) = R6(t - t) 

and 

E{w(t)v r (t) ) = E(v(t)x T (0)) = EfwCtJx 1 (0)} = 0. 

Q is a positive semi-definite diagnonal matrix and R is a scalar 
constant. The (4 x 1) matrix x(0) is the value of the state vector 
x at time zero. The estimator equations are then 

x = Fx + Gu + £H T R (y - Hx) , (3.17) 

where £ is the covariance of the error (x - x) which is determined 
by finding the steady state solution to the Riccati equation 

£ = FE + £F T + LQL T - £H T R _1 H£ . (3.18) 


It can be seen by dividing Eq. 3.18 by R 
state solution of £ depends solely on the matrix 


that the steady- 
Q/R. Thus, for 
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q 0 0 0 

0 q 0 0 

0 0 0 0 * 

0 0 0 0 

. — d 

the amount of work required to solve Eq. 3.18 is greatly simplified. 

Equation 3.18 was solved by numerical integration for a wide 
range of the parameter q and the results are shown in Fig. 11. In 

t —It 

this plot, [k^, K » K^, = [£H R ] , the gain matrix of the 

-4 -2 

filter. For driving noise with variance v = 10 sec and measurement 
-3 

error to = 10 radians, q would be 0.01. The solutions of Fig. 3-11 

are for parameter values of D = 1.5 sec 1 mean orbital rate and 
-3 -1 

n = 1.09 X 10 sec , corresponding to the 15 rev/day. Changing n 
-4 -1 

to 2.18 X 10 sec for 3 rev/day, changed the steady-state gains by 

less than 2%. Thus, there seems to be no need to use time-varying 
gains computed by continuous integration of the Riccati equations to 
account for the change in the orbital rate due to the elliptic orbit. 

Since the input to the Kalman estimator is not continuous in 

nature but comes as a sampled signal, the familiar Shannon sampling 

theorem requires that this signal have a sample rate which is at least 

twice as fast as the variation in the state which one is trying to 

produce. The rate terms a and a oscillate with a frequency D, 

x y 

therefore it is necessary to sample at least 21 / I times per satel- 

Z Z XX 

lite revolution. Thus, the horizon sensor heads must either be spun 

at least 21 /i times faster than the satellite spin rate or more 
zz' xx 

than 21 / I sensor pairs must be used on the satellite, 

zz' XX 

4. Magnetic Attitude Control of the Symmetric Spinning Vehicle . 

Given the ability to determine the three unknown states of the 
vehicle attitude, one can proceed to develop a control system making use 
of these states. Then, once a control law has been formulated, the 
implementation of the required control torques must be considered. 

The control problem studied here is that of a regulator. It is desired 
to keep the spin axis as close as possible to the orbit plane normal 
and maintain the spin speed of the satellite within some acceptable 
bounds . 


LQL of the form 


Q 

R 
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FIG. 3-11. KALMAN FILTER GAINS VS NOISE 
VARIANCE RATIO FOR SPINNING 
SATELLITE . 
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In this part, a new method of magnetically controlling the 
attitude of a spinning spacecraft is developed. It will be shown that 
this magnetic attitude control system is simple and requires very lit- 
tle power requirements for geodesy satellites in orbits between 15° and 
75° inclinations. This'" system, possibly supplemented with a passive 
nutation damper, thus becomes a strong candidate in any geodesy mission 
for which active attitude control i£ required. 

a. Magnetic implementation of the control law . 

The desired control torque will usually consist of components 
about the roll and yaw axes to correct wobble motion and pointing error 
and a component about the spin axis to correct spin speed error. To 
implement a torque magnetically, one makes use of the relationship 

T = in XB (3.19) 

where m is the magnetic dipole created by producing electromagnets 
or passing current through coils fixed to the satellite. The vector 
B is the local value of the earth's magnetic field and T is the 
resultant torque. 

Because B has an arbitrary direction (see Fig. 3-12), it is 
not always possible to solve Eq. 3,19 for m to produce a desired T, 
Thus, it is assumed that the spin component of desired torque is ignored 
except in cases where the spin speed has deviated so far from the 
nominal value that corrective action must be taken. 

For maximum efficiency, it is necessary to create a magnetic 
dipole in the spacecraft normal to the earth's magnetic field at any 
given instant. Thus, 


m • B = 0 . (3.20) 

By solving Eqs. 3.19 and 3.20 with the assumption that the component 
of T along the spin axis is zero, one obtains 

m z = < T y D B * ‘ T x D V/l B ! 2 

m x - <™z B x " T y D )/B z (3 - 

m y - ( m Z B y + T X D >/ B Z ’ 
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where T y and T Xd are the desired control torques and the subscripts 
x, y, z on all quantities denote measure numbers in the local (L) ref- 
erence frame. Equations 3. 20 is basic to magnetic attitude control and 
have been used in modified form in several mechanizations since appear- 
ing in Ref. 3-15. 

To mechanize a control system along the lines of Eq. 3.21 

requires the measurement of three components of the magnetic field. It 

also requires seven multiplications and two divisions. Although this 

is complex in itself, the worst problem arises from the requirements for 

division by and | b| F or a highly elliptic orbit (see Fig. 3-12), 

B alone varies over two orders of magnitude in size, 
z 


One method of control mechanization which greatly simplifies 


Eq. 

3,21 is to set m =0 
z 

and assume an average value of the 

spin 

component of magnetic field 

B . Then, for K = l/B , 

z ' z average 




m - -KT 

x yn 

(3.22) 



• 

EZ 

n 

>, 

s 


The 

resulting actual torque 

acting on the vehicle will then be 



T x 

= (KB z )T Xd 



T y 

= (KBz)Ty D 

(3.23) 


T z 

= " K(B y T y D + b * t x d > • 



The spin torque T should average to nearly zero over several orbits 

z 

under such a scheme. 


Pointing control can also be achieved with the z-coil alone. 
If one assumes that the magnetic field B is normal to the required 
torque T^ , then 

m = K 1 (B x T d ) 


where 


is some appropriate gain. 


Thus 


m z ~ K i (B x T y D B y T x d ) • 

Spin-coil control has the disadvantage that the resulting pointing con- 
trol torque is not usually in the exact direction of the desired torque 
but has the advantage of creating no undesired spin torques. 
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b. Spin control 


If the spin speed deviates from the nominal value, this can 
be detected by variations in the periodic signal coming from the horizon 
sensors. The desired control torque to correct such a deviation AW 

z 

would be 


or 



-K AW 
z z 



-K sgn(AW ) 
z z 


(3.24) 


when | AwJ exceeds some deadband value. There are two ways in which 
this spin control can be implemented. 


The first implementation of spin control follows from Eqs. 3.23. 
From these it is seen that a spin torque will exist because of the 
presence of pointing control. Thus, one might incorporate logic into 
the control system such that pointing control is only actuated when 
the resulting spin torque is in the desired direction. 


The second possible scheme is to apply current to the x and y 
coils in such a way that the resulting magnetic moment is normal to the 
component of the magnetic field in the X-Y plane. Thus 


m ‘ = KB sgn(AW ) 
x z y z 

(3.25) 

m = -K B sgn (Aw ) . 
y z x z 


c . Complete control logic 

The preceding two sections (a & b) , can be combined to yield a 

control system which provides both pointing control and spin-speed control. 

Two constants and are defined which represent boundaries of the 

deadbands associated with the spin speed. Let be the value such 

that when I Ad > C, , one should provide some sort of spin control. 

1 z' 1 

Let C z be the value for which, when | AW z | > C^, it is mandatory that 
spin control action be taken. Then suitable logic governing the currents 
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put into three orthogonal coils would be 


for: -C, ^ ^ C, 

1 Z 1 


(Mode 1) 


v 


m 

x 


-K T 


D 


m = K T 

y X] 

m = 0 

z 


(3.26) 


for: C n < Aw C (Mode 2) 

l z 2 


if sgn [ -K(B T + B T )] = 

y y^ x 


if sgn [ -K(B T + B" T )] = 

B L y y xx 

J D 


m =0 
x 

m = 0 

y 

m = K, (B T -BT) 

7 1 -V - TT IT V 

■'D ' “D 

for: C ^ | Aw | (Mode 3) 

a z 


sgn(-AW z >, use Mode 1; 

-sgn(-AW ) 
z 

(3.27) 


m = K B sgn (Aw ) 
x z y z 


m 

y 


= -K B sgn (AW ) 
z x z 


m 

z 


K. (B T 
1 x y. 


D 


- B T 

y x 


D 


(3.28) 


If the gains are set properly and care is exercised in design and 
construction of the vehicle, Mode 3 will not be required during normal 
operation. However, it should be provided for use during initial 
spin-up or for the case where spin speed changes but no pointing error 
accumulates . 
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This control mechanization requires the presence of three 
orthogonal coils about the three orthogonal control axes. The z-cdil 
(spin coil) could be provided with a ferromagnetic core so that weight 
and power requirements would be lowered. 

The control system must also have the ability to measure the 

magnetic-field components, B and B when in Modes 2 and 3. This is 

x y 

done with a two-axes fluxgate magnetometer [Ref, 3-16], a device with 
reasonable linearity. The magnetometer's sensitive axes must be 
mounted such that no interference is created from the magnetic field 
produced by the control coil. For the case where the magnetic compo- 
nents m and m are both being generated, this interference can 
x y 

best be prevented by a time-sharing procedure. Here, the coil currents 
are cut off momentarily once each cycle of spacecraft spin. At this 
time the magnetic-field measurements could be sampled and held during 
the controlled portion of the cycle. 

d. Skewed-coil magnetic control. 

It is possible to simplify to a great extent, the mechanization 
of the control system suggested in the previous section. This can be 
done by replacing the three coils by a single coil skewed at 45° to 
the spin axis as shown in Fig. 3-13. A further simplification is to 
use only one magnetometer with its sensitive axis along the node of 
the skewed coil and the' x-y plane as indicated. Such a magnetometer 
would not require time-shared measurements. For such a system, if the 
current in the coil is constant, an average magnetic dipole moment is 
generated along the spin axis. If a constant current has its direction 
switched every 180°, the average magnetic moment is in the x-y plane 
pointing in the direction 90° from the switch points. Thus, one has the 
ability to generate all three components of the desired magnetic dipole 
in averaged over a spin cycle of the satellite. A disadvantage to such 
a system is that it can only generate the average required magnetic 
components in the x and y directions over one cycle. Thus, it does 
not have the ability to generate control forces which fluctuate faster 
than spin rate such as wobble damping terms. Therefore, the skewed- 
coil must be supplemented with a nutation damper. At any rate, the 


- 79 - 




FIG. 3-13. GEOMETRY OF SKEWED COIL IN SPINNING SPACE- 
CRAFT. MAGNETOMETER SENSITIVE AXIS LIES 
IN PLANE OF COIL. 



packaging advantages of a system with a single coil and magnetometer 
seem to make the skewed coil mechanization worthy of consideration. 


e. Minimum power optimal pointing control . 

The normal mode of operation of the magnetic control system 
will be Mode 1. It is necessary to specify what the desired control 
torque components T x ^ and Ty^ are for this mode and what appropriate 
gains should be. A reasonable criteria for selecting the control torques 
to be applied is to choose T x ^ and Ty^ which yields the response 
desired and minimizes the power used to achieve this. The expressions 
for the idealized pointing control torques will now be presented 
below. 


Power consumption is proportional to current squared and the 
magnetic moment of a coil equals the product of the current, the coil 
area, and the number of turns of the coil. Thus, for a single coil 
one should try to minimize the control current squared. For two 
orthogonal coils, minimum power is obtained by minimizing the sum of 

~ n-s-j 1 THmo i mi 2 S'*" 1 O n prnhl otn 

LixL oijuux co uJ. liiC culiowi-o . y r ~ i- — 

considered is that of minimizing the performance index J subject to 
the system differential equations 3.4, where 


- f «, 

—OO X 


2 2 2 2 2 2 
a x + q i a y + q 2 0 + q 2 Cp + T D + T D )dt 

x y 


For a fourth-order system whose equations of motion have 
complex symmetry, the solution of this problem is well known (e.g., 
Ref. 3-17) and yields an optimal control of the form 



-K <t>-K.,<t>-K o e 
v pi p2 



-K 

v 


e + k 


P 2 



(3.29) 
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where the thi^ee coefficients of Eq. 3.29 are solutions to the equations 


2(Kpi> + (D + 2Dn + n + d + q^MK^) + 2(Dnq^ - q^K 
2 2 2 

“ (d q 2 + d n q 1 + q^) = 0 


Pi 


(3.30) 


K v = " d + ^d 2 + 2K pl + q^ 


K p2 = dn + y ( dn)2 + q 2 -<K pl > 2 - 2DnK pl . 

_ ^ 

The solutions to Eq. 3.30 for D = 1.5 rad/sec, n = 1.09 X 10 rad/sec, 
q = 0, with q treated as a parameter are shown in Fig. 3-14. 

12 _ 4_1 

The nutation damping coefficient, d, is varied from 10 to 10 . 

-4 

Changing the mean orbital rate, n, to 2.18 X 10 , corresponding 

to a 3 revs/day orbit had no markable change on the solutions. The 
choice of actual gains to be used for power minimization depends upon 
the coefficient of the nutation damper used and the response desired 
which is determined by choice of the parameter q . 


f. Stabilit y considerations 

In the previous section, gains from minimum power control 

were determined with the assumption that the resulting idealized torque 

would always be mechanized exactly. However, using the simplified 

control mechanization specified in Eqs. 3,26, 3.27, and 3.28, one 

would only be able to actuate the coils to produce control torques which 

would be optimum in some average time sense. For instance, if Eq. 3.26 

is being used and the gain K is chosen as the average value of l/B , 

z 

then the actual control torque produced varies considerably from the 
"optimum." It is known from optimal theory that the optimal gains, 
when mechanized exactly, will result in a stable system. But when 
the actually mechanized control torques are time varying, depending 
'upon the time history of B^, the question of stability remains 
open. The subsequent discussion is intended to be a summary of the 
stability investigation. A more detailed development may be found in 
Ref. 3-19. 
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FIG. 3-14. OPTIMUM CONTROL GAINS FOR POWER MINIMIZATION AS A FUNCTION 
OF DAMPING. THE PARAMETER q 2 IS CHOSEN FOR RESPONSE 
CONSIDERATIONS. 
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Consider the case of the skewed-coil controller. Let the 
damping terms of Eq. 3,29 be lumped in with the nutation damper so 
that the desired magnetic torques are 


r D = - K pl * - K p2 6 
X 


T D y = K p2 ^ K P 1 9 • 


(3.31) 


From Eqs . 3.23 and 3.26, the actual position control gains are 


T 

x 


- KCK pl ♦ + K P 2 0> V l > 


- K l< t >(K pl » + K p2 9); 


T 

y 


- K< - K p2 * + K pl 9) V‘> 


-K 1 <t»(- Kp2 ♦ + K pl 9) ; 


(3.32a) 


(3.32b) 


where K, (t) = KB (t). Substituting Eqs. 32 into Eq. 4 yields 
1 z 


r* 

a 

X 


-d 

-D 

-V t)K P i 

-V t)K P 2 


a 

X 

0 

a 

y 


D 

-d 

V t)K p2 

-V t)K pi 


°y 

• 

4> 

— 

1 

0 

0 

n 



• 

e 


0 

1 

-n 

0 


e 


(3.33) 


One must determine what range of K^(t) can be allowed and yet be 
assured that Eq. 3.33 remains stable. 


Referring to Fig. 3-14, it is seen that K 

K, (t )K by K and K, (t)K , by zero in Eq. 3.33. 
1 p4 1 pi 

proceed to determine what value of the constant gain 
the resulting system to be unstable. 


» K , , so replace 
Pi 

Then, one can 
K would cause 


For the values of d and n, typical of the satellite considered 
here, application of Routh's criteria produces the approximate condition 
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Dd > K x (t) . 


(3.34) 


Equation 3.34 gives an approximation to the upper limit on K^(t). 

As a more rigorous approach, a Lyapunov function has been constructed 
which demonstrates asymptotic stability for K^(t) less than an upper 
limit agreeing very well with Eq. 3.34. The system is stable for 
k^(t) < 5.398, for example, using the Lyapunov approach with parametric 
values, 


D 

= 1.5 


-2 

d 

= 1.26 X 10 


-3 

n 

= 1.09 x 10 


-5 

K pl 

= 2.85 x 10 

-3 

K „ 

= 3.5 x 10 

p2 



Again, the time variation of orbital rate has only a slight 
effect on the results. 

A method of determining the limits on the allowable position 
gain has now been established. Suppose one sets the gains of the 
control system so the system is optimal for the average magnetic field. 
Then the previous technique determines whether or not stability is 
maintained over the entire range of magnetic field variations for that 
particular orbit. The average magnetic field for an orbit can be found 
from plots such as Fig. 3-15 taken from Ref. 3-18 by using the semi- 
major axis as the average radius. 

One must also investigate the stability of Mode 2 control, 
i.e., the situation when only the spin component of the magnetic moment 
is used for pointing control. Consider again the case where the magnetic 
control is for position control only. By requiring that the average 
magnitudes of q^, a , 9, and $ be always decreasing, one can obtain the 
condition 

K (B 2 + B 2 ) 

d> ^5 (3 - 35) 

2 2 

and K =0 for all (B + B ) of the orbit. The spin magnetic 
pi x y 

moment during Mode 2 and 3 control should then be 
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k iV [ * b x + 0 v 


(3.36) 


m 


/ 

where is chosen so that Eq. 3.35 is satisfied. 


g. Simulation results . 

Throughout this study, analog and digital simulations of the 
system have been made to verify the analytical results. A brief summary 
of some of the findings of these simulations will be mentioned here. 

In deciding to use a state estimator as part of the control 
system, it was necessary to compare the system response with that which 
would occur when no estimator was used. Without an estimator the 
magnetic position control could be 



- k p e 



(3.37) 


where Q is the roll error signal coming directly from the horizon 
sensors. This control would, of course, be supplemented with a 
nutation damper. 

Assuming, the exact control torque is applied as specified by 
the control laws 3.37 and 3.31, the response times were compared for 
driving the states to zero from a variety of initial conditions. For 
the cases considered, the system with the estimator was found to be 
all the way from equally-as-f ast to twenty-times-as-f ast as the system 
with no estimator. 

This speed of response has two advantages. Because dis- 
turbance torques are going to move the satellite spin axis away 
from normal to the orbit plane, fast response means that the average 
deviation is less. Also, because the control torque depends 
upon the magnetic field available (which fluctuates considerably in 
elliptic orbits), it is desirable to apply control when the magnetic 
field strength is sufficiently high. However, without an estimator, 
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the detectable error (roll component) might only be present when the 
magnetic field strength is low. Thus, without an estimator (which 
serves as a memory of the error), it could take much longer to drive 
the system to the null position. 

Analysis and simulation showed no advantage of the estimator 
for offsetting a steady state yaw disturbance torque. But since actual 
disturbances (see Figs. 3-4 and 3-5) are fluctuating, the advantage 
is realized. The steady state error due to a constant yaw torque is 
approximately <t> ~ T /(Dn) where T is the normalized yaw torque. 

X X 

Other points studied by analog simulation were the effect 
of putting a sampled signal into the estimator and using the constant 
parameter n in the estimator equations rather than a time-variable 
f + q. Comparing a sampled input signal with the continuous estimate 

A 

of the roll error Q causes no problem if damping is done solely by 

nutation damper. But such a mechanization causes the estimates of the 

rate terms a and a. to have a phase lag behind the actual rates, 
x y 

The lag angle depends upon the speed of sampling. Figures 3-16a, b, 
show the response of the rates a x and a ^ when sampling rate is four 
times per second. It can be seen that damping is poor and finally 
nonexistant when o? x goes to zero. This problem can be rectified 
either by sampling Q before comparing with the input Q or by 

A 

putting a lead angle in the rate damping control signal. If 0 ! x 

and a lag OL and a. by an- angle A. then the damping control should 
y x y 

be 

T = -K (cos A a - sin A a ) 

D v x y 

x 

(3.38) 

T^ = -K (sin Aa + cos A OL ) . 

D v x y 

y 

Figures 3-16c, d show a and a using Eq. 3.38 with A = 60°. 

An uneven response was the effect of using a constant n in 
the estimator instead of a variable f + g (which is available from 
the horizon sensors) to account for orbital rate. The trajectories 
followed by the satellite spin axis were highly dependent upon initial 
conditions in such cases and often resulted in limit cycles about some 
point near the origin. No evidence of (instability) was seen from this 
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FIG. 3-16a and b; SAMPLING EFFECT ON ACTIVE RATE 

DAMPING. SAMPLE RATE: 4 TIMES/ 
ROTATION. 
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FIG. 


3-16c and d: SAMPLING EFFECT ON ACTIVE RATE 

DAMPING. SAMPLE RATE: 4 TIMES/ 

ROTATION. 
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condition, however. 

A final point which required investigation was the performance 
of the system during Mode 3 spin control. It can be shown that to 
keep the pointing error as small as possible due to the spin control 
torque, this torque should be applied over a 180° segment of the 
orbit. For greatest efficiency, this segment should be centered about 
the perigee point. To check the spin correction performance of this 
control and also observe the effect on the pointing error, several 
digital runs were made with the actual magnetic field and realistic 
disturbance torques included. 

Figure 3.17 is an example of the response of spin speed, roll, 
and yaw errors during Mode 3 control of a 15 revs/day orbit. For 
this example, the magnetic moment components are 


5 

m = -10 B sgn(,AU> ) 

x y z 


5 

m = 10 B sgn (/\£0 ) 

y x z 


0 

m = 3.5 x 10 (d>B + 0B ) 

z x y 


2 2 
measured in Amp - m . The magnetic field components are in Webers/m . 

The parameters d and -D were 0.0126 sec ^ and 1.5 sec \ These runs 

were based on boundary constants between Modes 1, 2, and 3 of = 

0.01 rad/sec and C = 0.02 rad/sec with nominal spin speed of one rad/sec. 

6 


B. TRANSLATIONAL CONTROL MECHANIZATION 
FOR SPINNING SATELLITE 


The function of the translational control system is to keep the 
satellite centered about the proof mass in the presence of disturbing 
forces. This is accomplished by actuation of gas jets based on an 
error signal derived from the position of the proof mass with respect 
to a null point fixed in the satellite. As previously pointed out, 
spinning the satellite attenuates the problem of placing this null 
point where there is zero mass attraction of the satellite on the proof 
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mass. Hence, spinning makes possible improved performance; however, 
it modifies the mechanization of the translational control. Assuming 
a linear controller, Lange [Ref. 3-20]showed that a simple correction 
term can be added to the controller for a non-spinning vehicle to 
stabilize a spinning vehicle. A practical mechanization of this has 
been simulated on analog and digital computers and built in our labora- 
tory simulator. Under controlled conditions, this mechanization works 
very well; however, stringent requirements on the electronic tolerances 
and the mass-center/null point alignment are necessary to prevent a 
certain phenomenon which we call "trapping". Accordingly, control 
mechanizations have been investigated which will alleviate the trapping 
problem in the simplest way possible. 

The Trapping Phenomenon 

The existing control mechanization for the two-degree-of-f reedom 
laboratory simulator consists of pulse-width pulse-frequency modulators 
with deadbands and lead compensation in each axis and the to x p cross- 
coupling terms. The deadbands along each body-fixed axis create a 
square aeaaspace in a plane. In Fig. 3-18 below, lines of control 
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direction are sketched in the first quadrant around the deadspace. 

Note that in the vicinity of the deadspace corner, the force is in the 
wrong direction except on the 45° line (it should point towards 
the origin). This fact, coupled with small errors in deadspace square- 
ness, cross-coupling terms, or mass-center alignment gives rise to 
control-force directions sufficiently erroneous to cause the proof 
mass to become trapped near a deadspace corner. Because the satel- 
lite is spinning, the mass center (origin in Fig. 3-18) is describing 
circles in inertial space. This state is stable and the control force 
continues with just the right direction and magnitude to balance the 
centrifugal force of the mass center's circulaz* motion, thus wasting 
propellant . 

The point where the center of the proof mass becomes trapped is 
readily calculated by assuming the control force is linear after 
crossing the deadspace threshold. The control law given by Lange 
[Ref. 3-20] is modified by the square deadspace yielding an expression for 
the control error signal valid in the shaded region of Fig. 3-18 as 
follows: 

e = x + y (x - toy) - d (3.39) 

X X 

e y = y + r(y + wx ) - d y (3.40) 

= rate gain/position gain 
= vehicle spin rate 
= deadband parameters . 

The approach taken is to assume the vehicle is in the trapped state 

(observed experimentally on the laboratory simulator and duplicated on 

the analog computer) and to look for points which will maintain it. 

The observed limit cycle consisted of constant values of e , e , 

x y 

x, and y with x = y = 0. Constant x and y imply the center of 
mass is traversing circles in inertial space if it lies anywhere other 
than at those values of x and y. Maintaining this state requires 
control action in the direction and magnitude to provide the centrifugal 


where 


7 

to 

(d , d ) 

x y 
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acceleration of this circular motion. 


Assuming the qenter of mass (c.m.) is located at (x , y ) , we are 

e e 

looking for a point (x,y) where 

e y - y 

V 6 

1 ) — = (direction requirement) 

e xx 

x e 


2) K 


/ 


e x + e y = \/(x - x )g + (y - y e > 2 


( 0 ‘ 


m 


(magnitude 

requirement) 


K = control force gain, 

P 

co = actual vehicle rotation rate, 

a 


m = vehicle mass, 
v 

The direction requirement, coupled with equations 3.39 and 3.40 with 
x = y = 0 , yield a circular locus of points with center at 


e = e 
x xc 


e - e 
y yc 


d 

y 


d 

x 


- 7 cod^ -(7 co +l)y g 
2u )y 

+ 7 cod -(7 2 cj 2 +1)x 
_ y e 

2 «y 


(3.41) 

(3.42) 


and radius = /e^ + e 

V xc 


2 

yc 


(3.43) 


Similar expressions can be derived for the other three possible 
quadrants. 


Trapping will be maintained at some point on the circle if there is 
a portion of the circle in the assumed quadrant (e^ > 0 , e y > 0 ), since 
values of K are typically large enough to satisfy the magnitude re- 

p 

quirement. A portion of the circular locus appears in the assumed 
quadrant when 


e > 0 

xc 

or (3.44) 

2 2 

d - yiod -(7 co + l)y >0 
y ' x e 

In the laboratory simulator, 7 C 0 = 1 and with x g = y g = 0 and 

d = d (square deadspace) there is no solution, However, d > d 
x y y x 

or y g < 0 , x g = 0 yield solutions. 
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The influence of the term 7(0 is interesting: 

710 < 1: there always exists a trapping solution in some quadrant 

even with a square deadspace and no center of mass error, 

yw > 1 : the requirement for perfect squareness diminishes. How- 

ever, the affect of c,m. error is amplified by the 
factor (72^2 +l)/ 2 yCJ 

It is also interesting to calculate the c.m. displacement (y ) at 

e 

which the trapping locus just appears: assume a square deadspace, i.e. , 


d = d = d 
x y 


then the locus just appears when 


or 


2 2 

yiod - (7 to + l)y 


= 0 



1 - 7W 


2 2 

7 w + 1 


( 3 . 45 ) 


Differentiating this expression with respect to yu) and setting equal to 
zero yields a maximum permissible excursion (y g ) with a perfect square 
deadspace 

* 

y 

= 0.206 


for 


710 = 2.4 . 


With a circular deadspace, the solution is no longer quadrant 
dependent, and the locus of possible solutions is again a circle 



( 3 . 46 ) 
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The, requirement for trapping is that the above circle diameter exceed 

the deadspace radius. Hence, the ratio of allowable center of mass 

excursion (r • ) to deadspace radius r ,, is given as 
cm db 


cm 


db 


ti iy 

~ 2~2 " 
CO y +1 


(3.47) 


Again, the maximum permissible center of mass excursion (T ) is readily 

cm 

calculated 


at 



(3.48) 


Analysis of an octagonal deadspace follows along similar lines as 
the round and square deadspaces and yields 


at 


* 

d 

cm 

___ 


0.45 


yCO = 2. 


(3.49) 


2. Control Mechanizations Investigated 

Basic ideas considered’ for the elimination or attenuation of the 
trapping phenomena have centered around modification of the deadspace 
shape and introduction of state estimation techniques. Desirable 
properties of estimators include 

(a) noise attenuation, 

(b) ability to estimate c,m. position (x e , y e ) and control to it, 
rather than sensor null. 

(c) high bandwidth enabling simpler modulators. 

Susceptibility to trapping is affected by (1) center of mass dis- 
placement from the sensor null point (x and y e ), and (2) errors in 
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mechanization of the control law. Results of analysis of the first 
category (carried out in Section 1 above) are given in Table 3-1. 

The maximum values of the allowable excursion are tabulated along with 
the value of yto which produces the maximum. Trapping, however, is not 
the only consideration in selecting the best value of 7 to for a control 
system design. The value of 7 (rate gain/position gain) must be 
compatible with good control system performance in the presence of large 
disturbing forces and large initial conditions. The angular spin rate 
(to) is influenced by attitude control and mass-attraction averaging 
considerations. These factors may dictate choice of a 7(0 value 

v 

allowing less c.m. excursion than that tabulated. 

Mechanization errors arising from component tolerances and ampli- 
fier biases typically cause control signals to be off about 10 % of 
full scale or about 20% of the deadband. This translates roughly 
one to one into an equivalent c.m. excursion, hence, a system could be 
on the verge of trapping from mechanization errors alone if a square 
deadspace was used. 

Table 3-1 also indicates relative complexities of the systems 
without compensation or estimator. The primary factor making systems 
(T) and @ simple compared to other schemes is the square deadspace. 
Overall system complexity is a strong function of the estimator con- 
figuration which will be discussed next. 

3. State Estimation Techniques . 

If a dynamical system is thought to obey the set of system equa- 
tions 

x = Fx + Gu 


where 


x 



A 

£ 


u 


A 


state vector, 
constant matrices, 
known control vector, 
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an "estimator" or "observer" for this system is 


x = Fx + Gu + K(y - Hx) 


(3.50) 


where 

x = estimated state vector 
A 

K = constant estimator gain matrix 

A 

y = measurement vector 
A 

H = constant matrix. 


When the measurement vector is also a state, the estimated state vector 
need not include these states and Eq. 3.50 is somewhat modified. This 
is sometimes referred to as a "reduced state observer" or "minimal-order 
observer". Recent references covering the theory of these state esti- 
mation techniques are Ref. 3-21 and 3-22. 


The normalized translational 
rotating body in body coordinates 


equations of planar 
may be written 


motion of a 
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— — 
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0 0 1 0 


X 


0 0 
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0 0 0 1 
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3 
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3 
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0 1 
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y 

• 

V — 


y e 


(3.51) 


where 

x,y 



f 


d 


x 


i 


X 



r: position of mass center with respect to inertial space, 

= time rate of change of position, 

= control accelerations, 

A . 

= spin rate, 

A 

= disturbance accelerations, 

= center of mass displacement. 
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If we assume 


d 2 

(a) disturbances and c.m. displacement negligible (f - to x = 0, 

etc.) / x e 

(b) x,y are measured directly, s 

then, one of the most simple examples of an estimator is example (jj) 
in Table 3-2. 



If f and c.m. displacement are not negligible, errors will be present 

A A 

in the estimate of the states (x and v). If these errors are not accept- 
able or if we desire knowledge of the unknown quantities, then the 
state vector is augmented with these errors as constants to be estimated. 
The estimator equations then become (example (e) in Table 3-2) 



(3.53) 


A reduced state estimator with a revised state formulation can also be 
mechanized for this system. If we define 
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Description 

Noise 

Rejection 

lead compensa- 
tion 

bad 

1(b) 

4-state (x, v) 
estimator 

best 

© 

2-state (v) 
estimator 
(reduced State) 

good 

© 

4-state (v, c.m.) 
estimator 
(reduced state) 

good 

© 

6-state (x, v, 
c.m.) estimator 

best 


C.M. 

Tolerance 

Requirements 


Number of 
Operational 
Amplifier 


best 

see Table 3-1 

2 

good 

see Table 3-1 

6 

bad 

see Table 1 

4 

bad 

best 

8 

good 

best 

8 
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x - toy 


V 

x 

V 

y 


4 

4 


y + tox ( 


inertial velocity in body coordinates, 


then the equations of motion are 


— — 


— 1 





- ~ 

• 

X 


0 co J 1 0 


X 


0 0 

V 


t r-4 

o 

0 

3 

1 


y 


0 0 

• 

V 

= 

o o » o co 


V 

+ 

1 0 

X 


i 

i 


X 



• 

V 


o 

3 

i_ 

o 

o 


V 


0 1 

y 


L i - 


y 


_ 



f 


y 


(3.54) 


The system matrix of Eq. 3.54 is partitioned as shown in line with the 
states being directly observed. This defines the following submatrices: 


11 


F 


12 



to') 

oj 

O' 

1 _ 



to" 

0 * 


G 


2 



0" 

1 


L 




K 


12 


K 


22 


estimator gain 
matrix 


and the reduced state estimator in block diagram form coupled with a 
simple pulse modulator is as shown in Fig. 3-19. 


This estiraator/modulator configuration is the rotating vehicle 
equivalent of the derived rate controller in common use today [Ref. 
3-23]. 


Other estimator forms considered include variations of the state def- 
initions as well as the number of states being estimated. In addition, 
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FIG. 3-19. ROTATING VEHICLE DERIVED RATE CONTROLLER 


various mechanizations are possible for a given estimator form. For 
example, the subtraction of (x - x) and (y - y) in full state esti- 
mators (e.g. Eqs., 3.53 and 3.54) can be performed in a separate amplifier 
or in each integrator of the system. A vector block diagram illustrating 
this mechanization difference is given in Fig. 3-20, a and b. 

Results of the analysis of the different forms and mechanizations 
can be summarized as follows: 

(1) State redefinition so that the error signal is a state: 

(a) reduces errors due to component tolerance 

(b) reduces number of amplifiers required; 

(2) Mechanizing extra amplifiers for computing the difference 
between measured quantities and estimated quantities in 
full state estimators substantially reduces errors due to 
component tolerances ; 

(3) Reduced state estimators are inherently sensitive to component 
tolerances because the step suggested in (2) above is not 
possible. 

Table 3-2 summarizes some of the characteristics of various 
estimator forms and lead compensation. All the estimator forms share a 
common characteristic response to a control pulse. The control signal is 
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FIG. 3-20a: ESTIMATOR WITH SUMMATION AMPLIFIER 



FIG. 3-20b; ESTIMATOR WITHOUT SUMMATION AMPLIFIER 
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fed into a velocity state integrator essentially duplicating the response 
of the actual dynamical system. Thus, one could say the bandwidth of the 
response to a control pulse is infinite. This is fundamentally different 
from lead compensation which relies on sensing a change in velocity from 
observed position data. This high bandwidth to control pulses allows 
use of simpler modulators with estimators. 

4. Conclusions 

Utilization of a square deadspace appears marginal even with per- 
fect c.m. location, eliminating schemes with square deadspace (l) 
and ( 4 ) in Table 3-1, 

Two desirable characteristics of the estimators are (1) noise 
rejection, and (2) ability to estimate c.m. location and control to 
it, rather than to sensor null, thus relaxing the requirement of physical 
control of the c.m. location. 

If an estimator is required, systems utilizing the high bandwidth 
to simplify the modulator appear to be the logical choice (schemes 
( 5 ) and © over ( 3 ) in Table 3-1). Although lead compensation may be 
adequate from a noise and c.m. standpoint, its selection over a simple 
estimator form (such as (b) or (c) in Table 3-2) is not immediate. The 
additional complexity in the modulator required when using lead compensa- 
tion is balanced by the additional complexity of the estimator. This 
choice essentially reduces to pulse-width pulse-frequency or "rotating 
derived rate". Unlike l/s plants, where the derived rate is clearly 
the simpler controller, the rotating plant requires additional complexity 
in the estimator making the choice an even one. 

Mechanization of a system with circular deadbands, lead compensa- 
tion, and a p-w p-f modulator (scheme ( 2 ) in Table 3-1) on our laboratory 
simulator is nearing completion. This should allow verification of our 
analog and digital simulations of the system. In addition, a six-state 
estimator (Form (e) in Table 3-2) is being built for substitution in place 
of the lead compensation (thus producing scheme (3) in Table 3-1). This 
will provide laboratory verification of the estimator concept. 
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C. INTERACTION BETWEEN TRANSLATIONAL CONTROL AND 


GRAVITY STABILIZED ATTITUDE MOTION 
/ 

This section gives attention to the analysis of a phenomenon which 
can occur with a drag-free gravity-stabilized vehicle. This situation 
could arise, for example, if the currently planned GEOS-C geodesy 
satellite was to be made drag-free by the addition of an "add-on drag- 
free device". The phenomenon is the possibility of attitude instability 
as a result of coupling (or interaction) between the translational 
control system and the gravity stabilized attitude motion. 

If the pickoff null and satellite mass center are not coincident, 
the pickoff interprets attitude motion of the satellite about the mass 
center as a change in the position of the proof mass, and the control 
actuators are commanded accordingly to null out this position change. 

If the actuators can exert a moment about the satellite mass center, 
then there is the possibility that attitude motion itself can generate 
attitude torques in such a way that the system is unstable. 


1. Nonlinear Equations of Motion . 

Consider a rigid, gravity-stabilized Drag-Free Satellite to which 

there is applied a collection of translational control forces, £Fj, where 

the subscript j denotes the control actuator. If R denotes 

c-pm 

the position vector of the satellite composite mass center (c) with 
respect to the proof mass (pm) [Fig. 3-21], we can then write 



R 

c-pm 



+ ZF. 



D 


(3.55) 


In Eq. 3.55, F and F pm , denote respectively, the gravitational forces 
G ^ 

acting on the satellite body and on the proof mass. F^ is the collection 

of external disturbance forces on the satellite, and in and m denote 

t pm 

respectively the masses of the satellite and proof mass. 


For the relative acceleration between the satellite's mass center, 
c, and the proof mass, pm, due to the gravitational forces, one has 
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FIG. 3-21: GEOMETRY OF THE PICKOFF NULL, SATELLITE 

MASS CENTER, AND PROOF MASS. 
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where m denotes the mass of the earth, and co is the (constant) 
© o 

orbital angular velocity of c about the earth. 


f g 


m m 


pm 


-co 2 R 
o c-pm 


/N /S 

3R R 
pm pm 


(3.56) 


where 1 denotes the unit diadic and R is the unit vector in the 

_ pm 

direction of R ( R = R /r ^ . 

pm \ pm pm pm / 


To write the scalar equations of motion of the Satellite, let 

A A A A 

f 0) = ( 0 _ , 0 f 0 ) be a set of reference orthogonal unit vectors 
i z j 

r _ 

centered at |_Fig, 3-22] ; 0^ is the unit vector in the direction 

from the earth's mass center to c, 0 ^ is the unit vector in the 

— , /N 

direction of the orbital angular velocity w , and 0 is such that 

O 

A A A A 

0 “ ^3 x • (9) is rotating in space with angular velocity 

co = co 0 Let (N] = n g» n 3 ^ 136 a set of orthogonal unit vectors 

centered at c and fixed in the satellite. The orientation of (N) 
is such that it is coincident, at any instant of time, with the principal 
inertia axes of the satellite. 


Any orientation of the satellite with respect to the reference 

✓S A A 


frame (0) is obtained by aligning (N) with (0) and then performing 

A 

three successive right-hand Euler rotations as shown in Fig. 3-22: 0 
about the axis aligned with n , 0 about the new direction of the axis 


aligned with n , and 0 about the final direction of the axis aligned 
2 3 

with n . This set of rotations produces the so-called nonclassical 
3 

Euler angles. 


A A 

The relation between the unit vectors 0^ and n^ (i = 1, 2, 3) 
can be immediately obtained as 
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FIG. 3-22: GEOMETRY OF THE ORBITAL AND BODY-FIXED 

REFERENCE FRAMES. 
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By noticing that 
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F . 
J 


4 f 


n +f n +f n „ 
jl 1 j2 2 j 3 3 


(3.58) 


and by expressing the left hand side of Eq. 3.54 in the coordinate system 
(0), the three following scalar equations describing the translation motion of 

c with respect to proof mass in the absence of external disturbances, F , 

are obtained 


m 


m 


X-2w Y-3w X 
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If r (- r n +r n + r n ) denotes the vector from the j trans- 
J v Jl 1 j2 2 j3 3 

lational control actuator to c, the moment T. about c , due to F., 
is then 
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T = r. XF. = (r F -r F )n + (r F -r F )n + (r..F -r.F.,)n 

j j J j2 j3 j3 j2 1 J 3 Jl jl J3' 8 jl j2 j2 jl' 3 


= T n +T n„+T„n„ . 

jl 1 j2 2 j 3 3 


(3. 60) 


If I , I , and I are the three principal moments of inertia of 
12 3 

the satellite (excluding the proof mass, which never touches the satellite), 
application of the angular momentum principle gives the three following 
equations, describing the rotational motion of the vehicle about c 


L"l + ['WV’K = 3 VVV Sln0 2 COS0 2 Sine 3 + E T ,1 

I 2 "2 + LWV’K = 3 “o(Ir I 3 )Sln0 2 COS0 2 COS8 3 + E T j2 
r 3"3 + (I 2' I l ) “l“2 = 3 “o' I r I 2 )COs20 2 Sln0 3 COS0 3 + E T j3 


(3.61) 


In Eq. 3.61, to , to , and to are the three measure numbers of the 
12 3 

components of the angular velocity of the satellite (with respect to 
inertial space), expressed in (N), and h is the constant angular 
momentum (with respect to the satellite body) of a yaw coupling rotor 

— A 

connected to the satellite in such a way that h = hn^ 

The expressions for an in terms of the attitude angles 0^ 

(i = 1, 2, 3) are 


(0 = (0 + CO sin0 )sin0 + (0 cos0 - co cos0 sin0 )cos0 

1 2 O 1 3 JL 4 O 1 ^ 3 

co = (0 + CO sin0 ) cos0 - (0,cos0^ - co cos0 sin0 )sin0 (3.62) 

2 2 o 1 31 2 o 1 2 3 

co = 0 „ + 0 , sin 0 „ + co cos 9 .cos 9 . 

3 3 1 2 o 1 2 

2. \leasured Quantities, Control Forces, and Torques . 

The quantities sensed by the proof mass position sensors are the 
three measure numbers of the components of the vector from the "pickoff 
null" (pk) to the proof mass (pm). Let g^, g^, and g^ denote the 
three measured quantities (i.e., ^p^-pm = °1 + g 2 n 2 + 63 *^ 3 ) and 
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e , e„, and e the three measure numbers of the components of R 
1 ^ 3 


expressed in (N) . Then 


pk-c 


R , = R , + R 

pk-pm pk-c c-pm 


and it follows that 


■) 

<K) 

>-■ 

) 


" e l 

_T 

"X" 

g 2 


e i 

+ c 

Y 

i 

crq 

w 

1 


_ e 3_ 


Z 


(3.63) 


where the superscript T denotes matrix transposition. The quantities 
e , e 0 , and e„ can be thought of as bias in the pickoff signal due 
to non-coincidence of the pickoff null and the vehicle mass center. 


The three measure numbers f . (i = 1, 2, 3) of the control force 

Ji 

applied tc the satellite by the control actuator are assumed, 

for purposes here, to be linear functions of the sensed quantities and 
their rates. 


Ef.. = -(k g. + k g.) (3.64) 

Ji P i v 1 

where k^ and k^ are "position" and "rate” gain constants respectively. 

For the measure numbers of the control torque, ST., a linear 
relationship is assumed of the form 


ET , 


Ef 

jl 


jl 

ET „ 

= L 

Ef 

j2 


J2 

ET 


Ef o 

J3_ 


J3 


(3.65) 


where L is a constant 3X3 matrix of effective force moment arms 
about the satellite mass center, c, and is generally non-zero due to 
errors in alignment and physical constraints on placement of actuators. 
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3, Equilibrium Orientations and Positions 

To determine equilibrium solutions of Eq. 3.59 and 3.61, assume 

0. = constant =0 (i = 1,2,3) 

l ie 

X = constant ^ X 

e 

(3.66) 

Y = constant ^ Y 

_ e 

Z = constant ^ Z . 

e 

By substituting Eq. 3.66 into Eq. 3.59 and 3.61, and by assuming 
that all quantities in Eq. 3.66 are "small" (so that products of small 
quantities in Eq. 3.59 and Eq. 3,61 may be dropped), yields the equili- 
brium solution 



Also assumed in this solution is that Ik I » 3co m, , (i.e. the 

1 p o t 

translational control natural frequency is very much higher than orbital 
frequency) a condition actually satisfied in practice. 

Equations 3.67 simply state that in equilibrium, without external 
disturbing forces, the proof mass is coincident with the pickoff null 
and the body principal axes and aligned with the orbital reference 
frame. 
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4. 


Stability 


To a first approximation, the measured quantities g (i = 1,2,3) 
can be expressed, from Eq. 3.63, as 



(3.68) 


where the subscript s denotes a perturbation from the equilibrium 
solution. 


Assuming that there are two translational control actuators for 
each axis of the satellite, mounted in opposition to each other, the 
matrix L can be written as 



(3.69) 


Linearizing Eqs. 3.59 and 3.61 about the equilibrium solution given 

in Eq. 3.67, and writing the result in the Laplace transform domain 

(1 eaving out initial conditions) after taking into consideration that 

I k I » 3w 2 m , the following equations of motion are obtained 
1 p‘ o t 
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Equations 3.71 and 3.72 represent the full six-degree-of-f reedom 

system for small motion about the equilibrium solutions given in Eqs. 

3.67. It is obvious from the equations that the parameters e e , e 

X Z o 

(pickof f null offsets from the satellite mass center) and £ , £ £ 

X Z *5 

(effective moment arms of the control actuators about the satellite 
mass center) cause coupling between the equations of translation and 
attitude. The stability of these equations as a function of these 
parameters is the primary interest here. 

To simplify the analysis initially, assume that e and £ are zero 

o o 

and investigate stability of the motion in the plane of the orbit. The 
equations for this may be written 


[V 2 + 2? g “g I 3 S+3 “! <I 2- I 1 ) + V Sk v K Vl + V2 , ] <5 3s + 

- K + s \) (v* - vJ = ° ; 


(3.73) 


/ 2 
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sk \ 
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p j 
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l o t j 
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| 3S 
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sk 
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1 Y + 
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< s 

IV P 
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3s 


The characteristic equation of this set of equations can be approxi- 
mated very closely by .the equation 


K(p + k)p 


1 + 


/ 2 r 2\ ir 2\ 

p(p + _p + r j - G(—P + r j 


= 0 


(3.74) 


( p2 + Tr p + r2 ) ( 


■p + r | (p 2 + 2S C P + 1 


„ „ A Wi + W 

where K = 


Vg 


A 


3(12 " II) 


1/2 


) 


A 

r = 


/ k /m 

V p' t 




k = — 2- 

k v W G 
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p = co G s 


G t 2^2. ^l e 2 ~ ^2 e 1^ 

W G ^1 € 1 + ^2 S 2^ 


The parameter £ has been included in Eqs. 3.73 and 3.74 to 
0 

account for some rate-dependent damping of the gravity stabilized attitude 
librations. This damping was not included explicitly in the original 
equations because the form in which it enters the equations is so depend- 
ent on mechanization. 

Viewing the characteristic Eq. 3.14 as a function of two parameters, 
K, and G, it is possible to sketch the loci of the roots as a function of 
K while holding G constant. There are, of course, six poles for the 
loci, four corresponding to the X and Y translational mode with normalized 
undamped natural frequency of r = \/ k p/ m t/ w G> and two corresponding to 
the gravity stabilized pitch attitude with normalized natural frequency 
of 1. There are five zeros for the loci, one on the negative real axis 
at -k, two very near the translational mode roots (the exact location 
depends on G), and two very near the origin on the real axis. 

Fig. 3-23 is a plot of the roots of the system calculated directly 
from Eqs. 3.73 (not the approximate characteristic equation, Eq. 3.74) 
for the parameter values indicated in the figure. Since realistic values 
of the frequency ratio r are around 100, it is necessary to show the 
loci with two different scalings, a macro-view to illustrate the behavior 
of the roots of the translational mode and a micro-view to observe the 
behavior of the gravity stabilized roots very near the origin. For the 
case illustrated (G > 0) it is seen that the loci cross the real axis 
relatively quickly for values of K > 0, but that for K < 0, the system 
is stable until relatively large negative values of K are reached. A 
similar loci plot for G < 0 would show that the zero on the real axis 
near +2 moves to the left-half plane and the loci again cross the 
imaginary axis for relatively small values of K, this time for K < 0. 


To obtain an analytical relationship among the four parameters 

«,€,£, and l and system stability, the Routhian array may be 
1 2 1 2 

constructed from the approximrte characteristic equation, Eq. 3,74, and 
conditions for imaginary axis crossover deduced from this. This yields 
the stability criterion, ' 
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FIG. 3-23. ROOT LOCI A3 A FUNCTION OF k FOR TRANSLATIONAL AND ATTITUDE 
MOTION IN THE PLANE OF THE ORBIT (G held constant). 


£ G^/si 3 ( X 2 " J l> 


(3.75) 


(£ e - JL e ) 
'' l 2 2 1/ 


for 


,/ '3 ' 

( v 2 - W 

V 3 ^ - v 

( Vi + v 2 > 


This relationship governs the imaginary axis crossover for the loci 

originating at the gravity stabilized poles and is the primary constraint 

on specification of the allowable values of pickoff null offset (e and e ) 

1 2 

and control actuator moment arms (£^ and £ ). From Eq. 3.75 using approxi- 
mate GEOS-C parameters 

<V 2 - Vi> •= °- 125 


(£ , fi 0 , e . i e „ in meters) for stability. This agrees very well with 
X £t JL A 

values obtained from the calculations of the roots of Eq. 3,73 used in 
plotting the loci of Fig. 3-23. 


Even though one must be careful not to violate Eq. 3.75 in the vehi- 
cle design, this requirement does not seem to be a difficult design con- 


straint , 


For example, with = + 1 cm, and 


6 = + 1 cm, we can safely 

£ 


have £ and £ = + 1 meter for = 0,05. 

12 — Ci 


The question of stability in roll-yaw and the effect of non-zero 
e and £ has yet to be fully investigated, but considerations as above 

O *5 

are expected to yield similar results : tolerances on pickoff null to 

mass center and control actuator moment arm specifications which are 
quite comfortable from the design point-of-view. 
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APPENDIX A 


The plots of this Appendix show acceleration of mean anomaly 
as a function of inclination of 450 km and orbital frequency, s, 
of 7 rev/day. The degree, & , ranges from 7 through 15. 
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APPENDIX B 


The figures presented in this Appendix are of acceleration of 
mean anomaly as a function of inclination for perigee altitude of 
450 km and s = m of 3 through 6 .and 8 through 12. The degree, l , 
is varied from 7 through 15. 
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APPENDIX C 


For this Appendix, acceleration of mean anomaly is shown for 
= 14, perigee altitude of 350 km. 
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DOCUMENTATION OF 

DRAG-FREE SATELLITE FUEL COMPUTATION PROGRAM 


PURPOSE 


DFSFUEL is a digital , program written in Fortran IV, level H, version 
II, which computes the translational control propellant required for 
sustaining a Drag-Free Satellite in a prescribed orbit for one year. In- 
cluded in the computation are the effects of 

1) atmospheric density 

2) control system limit cycle 

3 ) ti-ansient recovery from external perturbations 

4 ) solar radiation pressure 

5) vertical offset of mass center from position sensor null 

6) effect of upper atmospheric winds, 

PROCEDURE 

The heart of the program is the dynamic, 2-degree-of-f reedom atmospheric 
model of Jacchia [Ref. l] as modified by Keating and Prior [Ref. 2j for asym- 
etry of the diurnal bulge. For computational reasons, it has been found con- 
venient to use the 2-degree-of-f reedom polynomial fit of this model given by 
Sorenson [Ref. 3 ] rather than use the original equations of Jacchia which re- 
quire an integration. The polynomial agrees with the model to at least 3 
places so it is adequate for present purposes. The program steps the satel- 
lite around in the given orbit calculating density at each step and numerically 
integrating the drag to get the total drag impulse for the orbit. To this is 
added "worst-case" impulses for limit cycle, solar radiation pressure, transient 
recovery, and vertical bias in pickoff null. 

To account for the possibility that satellite surface forces need not be 
aligned with the control thrustors, a factor of V 3 is used in the impulse 
computation. A further factor of conservatism is provided for in a special 
factor called SAFTY by which the total impulse is multiplied, before computing 
the propellant mass. 

VARIABLES 

Listed only will be the variables which can be controlled by input. Most 
of them have default values assigned internally if they are not overridden by 
input specification. 
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INC 

INCO 

INC I 

INCF 

EP 

EPO 

EPI 

EPF 

HP 

HPO 

HP I 

HPF 

S 

50 

51 
SF 
ISP 

CD 

MASS 

AOM 

DB 

DVMIN 

VOFF 

DIST 

NUMDIST 

RSTIM 

OMG 

W 

CWIND 

BLGO 

FIO 

FIOB 

AP 


orbit inclination in degrees (no default) 
initial orbit inclination 


in cases where several 
different inclinations 
are to be examined. 


in cases where several 
different eccentricities 
are to be examined. 


in cases where several 
different perigee altitudes 
are to be examined. 


in cases where several 
orbital frequencies are 
to be examined. 


orbit inclination increment 
final orbit inclination 
orbit eccentricity (no default) 
initial eccentricity 
eccentricity increment 
final eccentricity 

perigee altitude in km (no default) 
initial perigee altitude 
perigee increment 
final perigee altitude 

orbital frequency in revs/day (no default) 

initial orbital frequency 

orbital frequency increment 

final orbital frequency 

effective I sp of the propellant in sec (default value = 30 
for nitrogen plus tanks) . 

satellite drag coefficient (default value = 2.2) 
satellite mass in kg (default value = 86 . 5 ) 

satellite area to mass ratio in m2/kg (default value = 1.2^”^) 

controller deadband in m (default value = 1.5\ - 3) 

thrustor minimum velocity change in m/sec (default = 4.4\ - ^) 

vertical offset from pickoff null to mass center in m (default 
value = 1 . 

magnitude of impulsive transient disturbances in m (default = 2^ ~ > ) 

number of transient disturbances (default value = 1. o\ 4 ) 

time to damp transient disturbance in sec (default value = 30 .) 

satellite orbit right ascension in degrees (default value = -4.6) 

orbit argument of perigee in deg (default value = 270 .) 

ratio of rotational rate of upper atmosphere to earth's rotational 
rate (default value = 1.3) 

angle between Greenwich and Equinox at time 0 in deg (default 
value = 0) 

-22 2 

index of solar decimetric flux in units of 10 watts/m /cps band- 
width (default value = 175 ) 

monthly average of F10 (default value = 175) 

index of magnetic activity in units of 2 ^ ^ gauss (default value 

= 10 .) 


- 173 - 



Appendix D 


DAY day of the year (default value = 170.) 

DF increment in true anomaly to be used in impulse computation in deg 
(default = 5 ) 

SAFTY factor of conservatism used in obtaining total impulse (default 
value = 1.2) 

PRINTlThese variables were included to give the option of plotting. 

PLOT (However, plot routines are not yet included so these variables 
should pot be changed from their default values. 

PROGRAM INPUT 

In general, the order and quantity of input data is arbitrary for any 
combination of the variables listed above. Any variable which is not explicitly 
given a value by input will be assigned its default value. To input a value 

I for a particular variable, the name of the variable is placed on the card (or 
card image) left justified to column 1 and the value is placed anywhere in 
columns 9 thru 20 as an E or F format number with explicit decimal point. 

Thus, there is always one variable name and value per card except when inputting 
a range of values for eccentricity, perigee altitude, orbital frequency, or 
inclination. In these cases, the name of the variable representing the initial 
value is placed on the card left justified to column 1, its corresponding value 
in E or F format in columns 9 thru 20, and the increment and final values in 
columns 29 thru 40, and 49 thru 60 respectively. 

The set of variables (HP, EP, s) is redundant in that any one ox them 
may be calculated from the other two, values must always be input for any two-- 
and only two--of them. This gives some flexibility in specifying the orbital 
parameters. (See samples on page 185 and 186 .) 


OUTPUT 

The output generated by DFSFUEL is generally self-explanatory. The 
only part which requires interpretation is the table labeled "default 
overrides." The order of the variables listed in this table is 


INC 

EP 

HP 

S 

INCO 

INC I 

INCF 

EPO 

EPI 

EPF 

HPO 

HP I 

HPF 

SO 

SI 

SF 

ISP 

CD 

MASS 

AOM 

DB 

DVMIN 

VOFF 

DIST 

NUMDIST 

RSTIM 

OMG 

W 

BLGO 

CWIND 

DF 

DAY 

FIO 

FIOB 

AP 

PRINT 

PLOT 

SAFTY 
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Whenever one of these variables has been specified as input, the value 
will be shown explicitly in the table, otherwise the default value was 
used and the symbol *** will be shown in the table. 

A sample output is included here to illustrate the result of running 
the sample input shown previously. 


REFERENCES 

D-l. Jacchia, L.G. , "static Diffusion Models of the Upper Atmosphere With 
Emperical Temperature Profiles," Smithsonian Contributions to Astro - 
physics , vol. 8, no. 9, Smithsonian Institution Astrophysical Ob- 
servatory, Washington, D.C., 1965. 

D^-2. Keating, C.M. , and Prior, E.J. , "Latitude and Seasonal Variations 

in Atmospheric Densities Obtained During Low Solar Activity by Means 
of the Inflatable Air Density Satellite," Space Research VII , 
vol. 2, North Holland Publishing Co., Amsterdam, 1967. 

D-3. Sorensen, J.A., Ph.D Dissertation, Dept, of Aeronautics and Astronautics, 
SUDAAR 380, to be published. 
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SAMPLE OUTPUT 


••»:•*»« II \*ir ^ r v>l| 1 KF iJ FUR URA.» H’fct SATelLME OPERATION FOR HNfc YEAR . 


OEFAULT VARIABLE OVERR IDES 


4S.OOOPOP 

P. I ooooo 









300.000000 

50.000000 

800.000000 




•0 ****••*«♦** 


100.000000 




6e«ec»ci»6«**e 




10.000000 




«»««««»«<••»$ 


200.00OO00 


/ 


t*$**C*4«*e t« 





SOLAR RADIATION 
235, 966 


FIXED WORST CASE IMPULSES (N.SECI 
L1HIT CYCLE 
IS* 263 


TRANSIENT RESPONSE 
66*667 


INCLINATION - 65*00 DEG 

ECCENTRICITY » 0.10000 

PERIGEE ALT (KNI REVS/DAY CM OFFSET IMP (N.SECI DRAG IMP (N.SECI TOTAL IMP (N.SECI MASS OF PROP ♦ TANAS (KGI 


300.000 


13.606 


61.732 

15216.930 

, 15596.550 

63.659 

350.000 


13.653 


60.366 

6361.236 

6739.668 

27.506 

600.000 


13.306 


59.037 

2817.167 

3196.099 

13.037 

650.000 


13.156 


57.768 

1335.326 

1710.968 

6.986 

* 00 . 000 t 


13.016 


56.696 

656.115 

1030.507 

6.206 

550.000 


12*676 


55.281 

3X6.130 

709.307 

2.695 

600.000 


SM-: 11.TS5 


56.100 

•It 7*666 

569.660 

2.263 



’ 1 K\ \\ . . . .■ 


... •. .. • 


.... - • . 
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6 r .c.ooo 

12.600 

52.952 

97.491 

466.330 

1.912 

700.000 

12.466 

51.837 

55.512 

425.245 

1.736 

750.000 

1?. 335 

50.752 

33.033 

401.681 

1.640 

HOO.OOO 

12. 206 

* 9. 698 

20.511 

388. 106 

1.584 

COMPILE TIME* 
tSTOP 

tSTOP 

2.58 SCC, EXECUTION TIME- 

6*97 SEC, OBJECT CODE- 

18096 

avrts, array area- 
521. 

521* 

544 BYTE$«UNUSEO- 230560 BYTES 


i 

c . 

c. 


177 



/ 


APPENDIX D 


(Cont) 


PROGRAM LISTING 
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/ 


1. //OFSORAG JOB IJ363, 616, 1,21, °A FLEW IMG* *WSGLEVEL® L 

2. //JOBL1B DU 0SNARE*SYS2,PR0GLI8»0ISP«(0LD»PASSI 

J- //KENT EXEC WATFOR 

6. //GOoSYSIN DO ® 

5. ANATFOR 

6. REAL*8 VARI301/* INC *,*EP »,*HP •«•$ *»MNC0 

To 1 • t • EPO '.‘NPO *,*SO ‘.MSP *,*00 •,‘HASS 

B. 2 •,«AOM •,•08 •,'DVNIN . 1 ■ 'VOFF S'DIST •.•NUKDIS 

9. 3 •♦•ftSTIR S'ORG *.»H A, 8 CL GO •••CWIRO *.*0? 

10. 6 • , • DAY «,*F10 •♦•PIGS *,*AP *,*N0 PRINT •♦* PLOT 

11. 5 •♦•SAFTY •/ , 

12. REAL»8 RL02,&L05,&L09,DSFIT1,DEFI.T2,DER.TJ 

13. REAL«S OgXP,OtOGiRLD12 

16. REAL«8 RLC3,RL06,RUK>,IU.07,RL08,RL01Q»TRP,RLO,LRO,RHl 

15. REAt.SB OMTOIPIS JSI/3C»1,010/ ' '<• 

16. REAL ISP. HASS, MtKOlS,i:;CO.IKC!fl»CF l KtlAO»lRRC,lHPYR, 

IT. UNPSOL,!fiC,INCD,l«PDG.B,L!S5G.fiAST0T,IMPY0T,lSPC«0 

IB. INTEGER PRINT, PI.OT.UPDN 

19. C 

20. C INITIALIZE SATELLITE PARAMETERS TO DEFAULT VALUES . . . 

21. C 


22. ISP • 30.0 

23. CD » 2.2 

26. HASS » 86.5 

25. AON . 1.2E-2 

26. C3 • 0.0015 

2T. DVMIN « 6.6E-6 

28. VOFF » 0.01 

29. OIST o 2.0E-3 

30. NUNDtS « 1.0E6 

31. RESTIM - 30.0 

32. C 

33. C IMSTULIIC serai::: CSS1T iCLAS rA£iT5TC« TC ?.*iS5J 

36. C 

35. ONG - *6.6 

36. w » 270. 

37. CWIND - 1.3 

38. BLGO • 0.0 

39. F10 - 175. 

60. OF • 5.0 

61. F10B- 175. 

62. DAY . 170. 

63. AP - 10. 

66. PRINT « 2 

65. PLOT • 1 

66. SAFTV • 1.2 

67. KOPT • 0 

60. C 

69. C START MAIN PROGRAM BY READING IN ANY DEFAULT OVERIDES 

50. C 

51. 1 READ I5,100,EN0«2 I DEFLTl.DUNl ,DEFLT2,0UN2 .0EFLT3.DUK3 

52. 100 FORMAT I3IA8,G12I) 

53. DO 3 J-1,30 

56. IF (0EFLT1 .EO. YARIJII GO TO 1 6, 3 ,6. 7, 8, 9, 10 , 1 1 , 12, 1 3, 16, 15 

55. 1 ,16, 17, 18, 19, 20, 21, 22, 23,26, 25, 26, 27, 28, 29, 30,3l,32,33),J 

56. 3 CONTINUE 

57. GO TO 1 

58. 6 INCO « DUM1 

59. INCF • OUR I *0.01 

60. INC! • 1. 
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ft l • 


•)uino**tn « INCO 

ft?. 


0.1 I.) I 

"J. 

5 

IPO - OUMl 

6A. 


fPF » OUMl « t'.Ol 

6'.. 


fpl • .1 

6ft. 


*<n»r • kupt ♦ 1 

ft*. 


'lUTfVlMWI . EPO 

ft*. 


00 f'» 1 

ft'*. 

ft 

MPv) « DO’M 

/o. 


HPr - UUMl ♦ I. 

u. 


HP! • 100. 

7?. 


KOPT • KUPT ♦ 2 

71. 


OUTOUMm « HPO 

7*. 


GO TO 1 

75. 

7 

59 - 0UH1 

7b. 


$F - OUMl ♦ .0.1 

77. 


SI - 1.0 

7a . 


KOPT • KOPT ♦ 3 

79. 


OUTnilMlM • SO 

80. 


GO T .9 l 

*»1 . 

8 

INCO « OUMl 

«?, 


INCI * OUM? 

83. 


INC F • UUM3 

8*. , 


OUTOUM (51 - 1NC0 

85. 


OUTDU u l 6) « INCI 

8b, 


0UTr»UM(71 ■ I NCF 

d?. 


GO f!) 1 

88. 

9 

PP.9 * OU-l 

89. 


EPl * OUM? 

90. 


EPF » 0UM3 

91 . 


KOPT * KUPT ♦ 1 

9?. 


OUTOUViHI - EPO 

93. 


UUU?UM(9) * fcPI 

95 . 


nutUUMUO) * EPF 

7i> • 


GO TO 1 

9b. 

10 

HPO a OUMl 

9 7. 


HPl a JIJM? 

9S. 


HPF * 0U M 3 

99, 


K0P1 * KOPT ♦ 2 

1 90. 


OUTPUMl] I ) - HPO 

101. 


OUTOU'M Wl • HP | 

10?. 


nuTouMt i3i « hpt 

101. 


GO TO t 

1 H. 

11 

SO = 9U M 1 

105. 


SI * DU*? 

l'Jf . 


$F « DU M3 

197. 


KOPT » KOPT * 3 

108. 


OUT DU* ( 1 ft 1 • SO 

10 »• 


0UTi)UM(l5l - SI 

1 !(.. 


OUTOUM(lb) « SF 

111. 


GO TO t 

11?. 

w 

I SP » DUM ] 

1 1 >. 


OU7DUMU7) • ISP 

115. 


GO 7 0 l 

IIS. 

13 

CO « UUM 

lift. 


DJTOUMI 18 1 * CO 

1 1 7. 


GO TO | 

Uh. 

1ft 

MASS • OUMl 

U9. 


0UTDUMI191 « MASS 

WO. 


GO TO 1 

I>1. 

15 

AO* • OUMl 
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IV, 

123 , 

I 

I2S, 
1 2*. 
127* 
123. 
I/'*. 
1 ) 0 . 

HI . 
132. 
HI* 
I »*. 
I IS. 
I -IS. 
137. 
lia. 
nv. 

HJ. 
HI. 
H?. 
H3. 
III. 
HS. 
IU. 
1-7. 
HS. 
I-*'. 
HI. 
HI. 
H2. 
IV*. 
H*. 
HS. 
HS. 
H7. 
I S'*. 
H'». 
I'»l. 
1M. 
1 * 7 . 
HI. 
!»•<.. 

I ' v. 
H'.. 
H'. 

Hi. 

HO. 

1 /« . 
in. 

I 7». 
t M. 

I 7*.*. 

I 7 S. 

I 7*.. 

• 

I »* . 

1 

I - . 
HI. 


* AOM 

GO TO l 
lo Utt « OUMI 

0»JIDU*12li • Ot) 

GO TO l 

17 OV**IN • DUHi 

CUTl»UM|??| » OVMIN 
GO TO I 

H V077 «* OUMl 

nuTUUM|23l * VOFF 
GO TO X 

\9 0 1 ST • OUMl 

UUTI>U*I24» * 01 ST 

GO TO I 

23 N010IS - otm 

OUT OU^l 2S ) * NUM01S 
GU TO I 

21 RSUH ■ OUMJ 

OUTUUSI?*) * RSTIM 
GO TO 1 

2? OMG « OUMl 

OUTnu^(27) = CMG 
GO TC | 

2 3 * • OUMI 

OJTDU*M 2 3 I = W 
GO TO 1 

2* BIGO * UJ“1 

PUTi)-JMI29) = HI GO 
GO TO 1 

?S CWlND * OUH| 

UUTUUMI30I = CW [NO 
GO TO l 
?6 1>F * DU*<1 

OUTOUMi 31 ) = OF 
GO TO l 

27 UAV * OUMl 

OUTPJMT32I = OAY 
GO in j 

2s fn * m>*i 

"U T HUM (331 * f 1 0 
GU TO J 

29 F 10B * OUMl • 

r*uT;*U M (3A I * HOB 

Go in i 

10 AS = OUMJ 

flUT i>U M t 3 S ) ■ Ap 
0 1 TO t 
M PR|*;T . 1 

**#. » * HkPJf 

; t t i 

*2 PI T * 2 

OU»l> \J W I 7*1 * PLOT 
C.l* TO 1 

13 jAf TY * OU*M 

MIUI'MI *K | * SAfTV 

0* 12 1 
2 S’ *• T [ *| j| 

IF (P-l*.T .r J. II v.« • tf 

• «MT r UtICII 

HI ««M*AT I • l • *//<*• ^ >’H L A *1 T HtJJlPfrU FOK DRAG F«tt SAlfU Mt (X»c* 
f . .*.i 1 z \ • ////////«dX • *1 f ► AUl T YAHlABLf UVlRRUtSM 
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IH1. 

1*4. 

1 

I 

I *1* 
1*8. 
1*9, 
l*J rt . 
191. 
1 Q 2. 

193. 

194. 

195. 
194. 
I Q 7, 
198. 
J99. 
200 . 
2QI. 
?«?. 

203. 

204. 

205. 
204. 
207. 
2 r *«*. 
209. 
2IC. 
211 . 
212. 

213. 

214, 
?1S. 
216 . 
217. 
2J6. 
21°. 
22 * . 
221 « 
222 . 
22 1 . 
229, 


/ > 7, 
22*»* 
2?'« , 
2 r . 

: *i , 
> 3’, 
.Ml, 

• H, 
/ 14. 
. 1*». 

* »/. 


24 , 


WAITE *6,1021 10UTDUMH) , 1«1#3S) 

102 f PRH AT «///«• • ,8(?X f f-13.6)/) 1 
47 CIA U SOL 
C 

C ASSIGN FIXED CONSTANTS ANO PERFORM DEGREES TO RADIANS CONVERSION 
C 

WE * 0. 72926- A 
TP * 0.0 

FO * 0.0 ' 

G*E « 3.9858E3 
RE « 6371. 

FSTOP * 360.0 
TO - 0.0 

UAYP * 4 . 

A S'J *= 1.496F8 

SUTNC * 23.45 

WSU « -77.59749 

KOPT » KOPT - 2 

IF <KOPT . EQ. 21 GO TO 80 

RPF * HPF ♦ RE 

RPO » HPO ♦ RE 

PP1 » HP! 

*0 CONTINUE 

SEC PY » 365. *24. *60. *60. 

KRAD *.8 

SORT * SQRT (3.0) 

OTR »1. 257.29578 
R TO = 57.29578 
1NC0 * OTR*lNCO 
I NC I * DTR* INC 1 
fNCF * OTR* INCF 
QMG *l)T R*O^G 
W * OTR* W 
Bl GO * OTP *BLGO 
FO « OTK*FO 
FSTOP = DTP4FST0P 
PI = 3. Kl 5927 
T P I = 2.*P1 

REVCCN a 40,*60.*24./TP* 

CO = C'JSICMGI 
SO = S1M0MC) 

PUT = PI/2. • 

TPOT = 3.*PUT 
>jb • ;)TK*jr 

c 

C PFIt-tMlNC $U*L TF*PERATUKt PARAMETERS 
r 

TF *3 = 4 l H • ♦ 3 . 6*F 1 Ofl ♦ l . 8 *{ F 1 0-F 1 03 1 ♦ I 0 . 3 7*0 . 1 4* S I Nl T PI *C D AY- 151 . 
I l/'o5.1 1 *F luc** SI N(7.*TPi*|UAT-59. 1/365. 1 
RK4 1 * * 45 . *1)1 P 
•*<\/ » 1 2 • ♦ OT R 

» ap -f <C?»-O.Cr.«A 3 J » 

c 

• CiAPJtt FULL IMPULSE fop LIMIT CYCLE. TRANSIENT RECOVERY, AND 
V'll' .>9.1, I4« StlMlV, CASt CONDITIONS) 

» 

UlU * *«,•,, I/.JVMIN 

r«*u * >.**ass* .‘VMr*'* r ,rcpY/ uuic 

I * *M '• * •1ST**. , .1Y*1IS*MAS5/»«FSI I.M 

|m*S 1L * Sjc l*9A5S*A'j J »4.5E-6*SECPY*RKAD 

1> I*'?*.' j, II , * Ti • / •» 
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?44. 


NRlU <6,1031 INPSOI* IMPIC. 1HPTR 

M*., 

103 

FORMAT l////44X,*Fl XtO WORST CASE INPULSES IN . S£C 1 1 // l 3X, • SOL M R A 

746. 


1 0 1 A I | 09 • ,27X, ‘LIMIT CYCLf » ,25 X, * TRANSIENT RESPONSE 1 //3(l5X»flO.J* 

?4 7, 


215X11 

743. 

6* 

continue 

?i*0. 

C 


?sr». 

C DETERMINATION OF SUN POSITION 

?M. 

C 


7*7. 


FPSU « l, 6725 92 E- 02 

25\. 


WSU *> otr*wsu 

254. 


SU1NC * JTR« SUING 

?ss. 


OF SU * 0. 9856091*0TR 

75b. 


S'JM » OFSUMOAY-OAVPI 

257. 


ESI* SUM 

75** 


0.1 41 1-1,5 

2*>9. 

41 

tSU ■ SUM ♦ EPSU*S1NUSUJ 

260. 


CES * COS USUI 

261. 


CF3 » ICES - EPSU1/CI. - EPSUPCESI 

2b7. 


SFS - ISIMESUI*SQR T ( l.-EP$U*EPSUII/I l.-fc PSU*CES> 

2b3. 


SwS » SlN(wSU) 

264. 


CWS * COSfWSU) 

2oS. 


SUS • SWS*CFS ♦ CWS*SFS 

244. 


CUS - C WS*C F S - SWS*$FS 

267. 


SUIM . ARSIN1SUS*SIN| SUINC) 1 

2*3. 


SUOL ° ARCUS(UJS/CUS( SUL*1 1 

»4«*. 


IF ISUS.LT.O.OI S'*U - TPI - SUBL 

27». 

C 


2>l. 

r initialize the farilv of orbits to be invest igateo. compute orbital 

2 72. 

c PA«A*F TfcRS.ANO THEN INCREMENT THROUGH THE FAMILY 

2 73. 

C 


274. 

64 

INC - INCO- INC 1 

774. 

52 

INC * INC ♦ INC I 

2 7b. 


IF (INC.GT.1NCF) GO TO 1000 

277. 

6S 

IF (PRINT .CO. 11 GO TO 69 

?t*. ' 


INCH * RTD*INC 

279. 


WP III (6,1041 INCO 

?ao. 

l rt 4 

FORMAT I/////7X,* INCt 1 NATION ■ «.F 7.2,* DEGM 

?dl . 

49 

CON II SUE 

?M7. 


GO Tf) (43,56,591 ,«'»PT 

73*5. 

43 

EP * FP3 - CPI 

734. 

44 

FP a fP ♦ fcPI 

>86, 


IF UP.GT.FPFI GU T" 57 

73b, 


IF (PRINT .fci. 1) GO 70 

2 i 7 . 


W«1H (6,1051 EP 

2 ■* 4 . 

l*'S 

f'Jw^AT (///7<, •cCCE. , IT?ICI Tr - •,FT,5) 

739. 


w° I T E 16,1071 

290. 

107 

format (/// 2X,*PfKlGFF alt (KM»«, 8X , c RF VS/O AY • • SX.»CH OFFSET IN 

29 1 . 


IP IN.SECJ*. 1 X , • jPAG IMP (N.SfcCI 0 * IX * * TOTAL l«P (N.SECI* i 1X#**4*S 

2 >7. 


? S OR PRnp ♦ T AN 4 S (*GI'//I 

7 *1. 

7n 

COOT iNUfc 

774. 


PP * RPQ - R P 1 

795. 

ss 

RP * KP ♦ PMl 

296. 


IF (RP.QT.PPF1 GO TO 44 

79/. 


S - R E VC UN* S'JRT ( G H L / ( p .P/( 1 . - EP||**3) 

793. 


GO 1? 6? 

799. 

66 

S * SO - SI 

300. 

47 

S * S ♦ SI 

301. 


ir is.3T.sm r,c ru 4,; 

302. 


IF (PRINT . F). 11 GO P ?1 

305. 


WRITE ( b , 1 06 1 S 

3"4. 

Pi 

r'JRMAl (/// 7X, »RF VIIUTIUNS PER OA Y ■ SF7.3J 
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3W, 


WAITE Ifc.lOSl 


$06, 

100 

FORBAT 1 III 2t, •ECCENTRICITY* , »*,• PER ALT IKKIS SX.'CX 

OFF St* , 

107. 


IT IBP IN.SECI*, l X # *ORAG IMP IW.SECI*. U,*t3IAl IBP CN.SEC 1 • . IX 

SOS* 


2*fe$tSS 0? PRO? 6 TANKS (KGI * // 1 



T1 

CO^Ttfv'UE 


HO* 


EP « EPO - EPI 


hi. 

58 

FP • EP ♦ EPI 


312. 


IF IEP.GT.EPFI GO TO 57 


313* 


RP ■ ( I .-EP )• (GRE*REVCQN**2/SP*2 I ••I1./3.I 


316. 


GO TO 62 


SIS. 

59 

S * SO - S! 


316. 

60 

S - S ♦ SI 


SIT. 


IF (S.GT.SF) GO TO 52 


318* 


IF I PRINT • EQ. U GO TO 72 


sn. 


WRITE (6,1061 S 


320* 


WRITE (6,109) 


S2U 

109 

FORMAT (/// 2X,«PER ALT CIWH, 8X ,• ECCENTRIC I IT' , 5X,*CH 

OFFSET II 

322. 


IP 1 N.SEC) * » IX, 'DRAG IMP 1 N. SEC )•, IX, ‘TOTAL IMP (N.SEC)', 

IX, *MAS: 

323. 


2 OF PROP ♦ TANKS (KG)*//) 


32%. 

72 

CONTINUE 


325. 


RP ■ RPO - API 


326. 

61 

RP » RP ♦ RP1 


327. 


IF (RP.GT.RPF) GO TO 60 


320. 


EP • W - RPP|SP62/IGHEPREVC0N*62nP6U.23.l 


329. 

62 

CONTINUE 


330. 


A * RP/ll. - EPI 


331. 


HP RP - RE 


332. 


QRF • SORT C OWE /A* 63) 


333. 


OOA « l./A 


336. 


PER « TPI/ORF 


335. 

63 

CONTINUE 


336. 


SINC - SIN(INC) 


337. 


CINC ® COS ( INC ) 


338. 


SOC INC • $0«C!WC 


339. 


COC INC * COPCINC 


360. 

C 



361. 

C THIS SECTION UPDATES THE ORBITAL POSITION OF THE SATELLITE 


362. 

C 

USING A CENTRAL FORCE FIELD MODEL. 


363. 

C 



366. 

60 

F2 - FO - OF 


365. 


T1 » TO 


366. 


IWPOG - 0. 


367. 

38 

F2 - F2 ♦ OF 


368. 


F » F2 


369. 


IF (F.GT.TPII F - F - TPI 


350. 


COSF o COSI F ) 


351 * 


CE - (EP ♦ COSFI/Cl. ♦EP*COSF) 


352. 


E » ARCOS(CE) 


353. 


IF (F.GT.PI 1 E * TPI - E 


3*6 1 


N - F - EPPSINIFI 


355. 


T * TP ♦ H/ORF 


356. 


IF CF2.GT.TPn T - T ♦ PER 


357. 


SINF ° SIN(FI 


358. 


R * A* ( 1 • * EPPCE) 


359. 


V - $QRT<GHEM2./e - 00A> I 


160. 


OUM * l./SQRT ( 1 • ♦? .* c P®COSF ♦ EPPfPI 


361. 


SG - EP«SINF«OUH 


362. 


CG • (1. ♦FP©CQSF|«OUH 


363. 


U * W ♦ F 


366. 


IF (U.GT.TPf )U • U - TPI 


165 • 


S»J - S INI U) 
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366. 

W. 
:«68. 
160. 
3 70. 

371. 

372. 

3 7J. 

374. 

375. 

3 76. 

377. 

378. 

379. 

380. 
381 . 

382. 

383. 

384. 

385. 

386. 

387. 

388. 

389. 

390. 
391 . 

392. 

393. 

394. 

395. 

396. 

397. 
39R . 

399. 

400. 

401. 
40? * 

403. 

404. 

405. 

406. 

407. 
403. 

409, 

410. 
411 . 

412. 

413. 

414. 
416 . 

416. 

417. 
414. 
419. 
470. 

421 . 

422. 

423. 

424. 

4 25 . 

4?6 . 


CU - COS(U) 

Ct - SOUTH* - I SU®$! KC 1® ®2 1 

1FCABSI1NC - POT 1 .IT. 0,0001 1 GO TO ?7 

IMU.LE.Pl > CU 10 TO 

BEAM * Qt*C *TPt - APCOSCCU/CLI 

GO TO 79 

78 BLA* * OMG ♦ APCOSCCU/CLI 
GO TO 79 

If IFCU.LE.POT) GO TO 34 
1FCU.GT.TP0TI GO TO 34 
OLA* * OMG 4 Pi 
GO TO 79 
34 BIAM « OMG 

79 IFCBLAM.GT.TPll 8L AM * RL AM - TPI 
VSG « V»SG 

VCG * V*CG 

MO • CCO*CU - SOCINC*SU)*VSG «tCO*SU ♦ SOCI NC*CCJI*VCG 

Y ID * CSO*CU ♦C0C!MC9$U)PVSG ~CSO*$U -COC INC.ACUI9VCG 

210 * S1NC*$U®VSG ♦ SIMC9CU9VCG 

VW * CW1ND*WE«R*CL 

SBL * S l N( BL AM 1 

CBL • COS ( BL AW I 

XRD » X ! 0 *VW*SBL 

YRO * Y | 0 - VW*CBL 

ZRO » 210 

St * SO* S INC 

l AM » ARSIN(SL) 

C 

C DETERMINE ATMOSPHERIC OENSITY VIA SORENSEN POLYNOMIAL 
C FIT 70 UACCHIA TWO PARAMETER MODEL 

C 

ETT - 0.5*A9S(LAM-SUtMI 
THT *0. 54AHSCLAM* SULM) 

HA *131 AM - SUBL 

TASU * HA -KK45 ♦ KK 1 2 *S I N ( HA* RK 45 > 

IF ( TASU.LT.-PI 1 TASU «TA$U*TPI 
IF CTASU.G7.P1I TASU « TASU-TP! 

1IFP «TCn*Cl. ♦0.28*S!N(ETT) * 0 .2 B* C COS l THT >-$ INC ETTI !• (COSC O.S* 
1 TASUI ) * * 2 • 5 1 
T IF =» T IFP ♦ AP* 

FX * CTIF-9 00. 1/1750. ♦I1.722E-041*CTIP~800. 1*42 1 
FS * 0.029l*EXPC-FX*F X/2. 1 
HE - R - RF/SQRTC 1 . ♦ ( 6 . 7385E- 03 ) *SL *SL 1 
TMP =TfF-(rfF-355.l*EXP(-FS*CHE-120.1) 

KRO = 1 

IF CHE.LT .120.1 GO TO 81 

IF CHfc.Lf .200.1 GO TO 42 

IF CHE. IT .300. 1 GO TO 43 

IF CHfc ,lT 1 GO TO 44 

IF CMF.LT.600.1 GO TO 45 

IF (Hf.LT .600.1 GO TO 46 

IF CHE.LT. 700. 1 GO TO 47 

IF CHF.LT.HOO.I GO TO 48 

IF IMF. LT. 900. I GO TO 49 

G'.l TO 50 

81 WHITE 16.110) Hf 

HO FliNMAT ( 5 A * • Ml I GH T * , ,F7.l* , !S OUT OF RANGE *7 
GO TO 73 

4? RLCJI2 * -11 .•AL0GC2.461I 
HLO = 120. 

K Ri^ - ? 
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*? T. 

KLO » RL012*1.2* 




43 RID? * -0» 30540*70* 02 ♦ THP* ( 0. 29 1 47450-02 *TMP 0 (-C 

). 10723230-05 



l ♦0.38342H9J-1 1* TMP) > 




IF (KK9.EQ.2I RHJ * P L02* i • 25 



4 31 . 

IF (KRO.EQ. 21 CO TO 51 



*32. 

HID « 200. 



433. 

KRO • 2 



434. 

RLO a Rl(l2 



43*. 

44 R L 03 « -0.3817 3060*02 * THP* ( 0 • 1 20 1 8970- 01 ♦ TMP 6 {-0 

1,63458190-05 


436. 

I *0. 1 l 50696D-0B*THP> 1 



437. 

IF ( K RO. E Q. 2 I RHI * PL03 



439. 

. IF (KRO.EQ. 21 CO TO 51 



4 39. 

HtO *=300. 



440. 

KRO * 2 



44|. 

RLO = RLU3 



44?. 

45 PL'34 * -0 .41947990*02* THP* 1 0. 1443 33 60- 01*1 MP4 ( -0. 

70408650-05 


443. 

1 *0. l2043800-08*TMPn 



444. 

IF IKR0.EQ.21 RHI -PL04 



44*. 

IF (KRO.EQ,? 1 CO TO 51 



446. 

HLO =400. 



447, 

KRO * 2 



448. 

RLO »RL04 



449. 

46 RL05 = -0.45682130*02 ♦ TMP*( 0. 1 7523010-01 *TMPP(-0 

• 8304390 0-05 


4*0. 

l *0. I 3963050-03* TMP ) 1 



4*1, 

IF (KRO.EQ. 21 RHI * PLUS 



4*?. 

IF (KRO.EQ. 2) GO *U 51 



4*3. 

HLO = 500. 



454. 

KRO * 2 



4*5. 

RLO = RL05 



456. 

47 RL06 * -0.49073530*0?*T«P*(Q, 203069QD-0 l *TMP* 1 -0. 

94092410-05 


457. 

1 *0. 15592820-03* TMPM 



450. 

IF (KRO.EQ. 21 RHI *RL06 



459. 

IF (KRO. t Q. 2 1 GO TO *t 



460. 

HLO =600. 



461. 

KRO = 2 



462. 

RLO * RLU6 



463. 

48 RLO 7 = -0.50092830*02 ♦ TMP* ( 0.204 71 030-01 *T*P*(- 

0.69300820-05 


464. 

1 *0. 14249719-08* IHP1 1 



465. 

IF (KRO.EQ. 2 1 ®Hl r RL 07 



466. 

IF (KRO.EQ. 21 GO TO 51 



467. 

HLO =700. 



468. 

KRO * 2 



46° . 

RLO * RL07 



470. 

49 RL08 * -0.49687380*02 ♦ TMP* ( 0 . 1 5525720-01 *TMP*(- 

0.56217010-05 


471. 

1 *0. 73357390-09* TMPM 



47?. 

IF (KR0.EQ.21 RHI *RU)A 



4 73, 

IF (KRO.EQ. 21 GO TO 51 



474. 

HLO =800. 



475 . 

KRO * 2 



476. 

RLO *R108 



477. 

50 RL99 = -0.45953530*02 ♦ TMP*( 0 .673299lD-02*TMP*(-0 

.26290550-06 


478. 

1 -0. 30583030-09*THP) 1 



479 . 

IF (KR 0.EQ.2) RHI -RL09 



480. 

IF (KRO.EQ. 21 GO TO 5L 



Ac 1 • 

HLO * 900. 



48?. 

RLO * RLCJ9 



483 . 

R1010* -0.42505970*02 *TMP*( -0. 50 2440 30-0 3*TM P*(0 

.37461300-05 


484. 

1 -0. 1018655D-08*TMP1 1 



485. 

RHI -RL010 



486. 

51 LRO » RLO *(HE-HLOl*( RHl-RLOl *0.01 



487. 

RHO *DEXPILR01*1.0k3 
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4*4. 

C 


409, 

C CiMPUft ATMOSPHERIC DRAG IMPULSE 

490. 

C 


4*U 


OfcCT « T - T\ 

492. 


VHAGS « <XP0»«2 ♦ VRD**2 * 2RD««2)*l.£6 

49i. 


IMPDG * 0.5*aH0»AC»<*MASSPVMAGS4CD«DELT4SQR3 ♦ JNPOG 

*94. 


7 1 * T 

495. 

39 

IF CF2.L T »F $TLP ) GO TO 30 

496, 


IMPUG * INPUG*S*165. 

497. 


VACC * 2.*IGMt/lA**3f|*VOFF 

4Vd. 


1 tJPCMO “ MASS*SECPY®VACC 

499. 


IWPTOT ■ IMPOG ♦ MPLC ♦ I MPT K MMPSOl ♦ IMPCNO 

590. 


MASror « lMpTl)T*SAFTY/( 9.041 SP» 

401. 


IF 1 PR ! NT .63, 1) GO TO 73 

50?. 


GO TU 1 74,75,761 ,KOPT 

503. 

74 

WRITE (fc 9 lil» HP # S, I«PCMO, IMPDGi I MPTOT ,HASTOT 

504. 

111 

FORMAT {• * ,6(5X,F9.3,6X1/) 

595. 


GO TO 73 

504. 

75 

WRITE (6,1111 EP,HP, IMPCMO, IHPOGt IMPTOT.HASTOT 

507. 


GO TO 73 

SOH. 

76 

WRITE (6,111) HP,fcP, MPCMO, IMPOG, IMPTOT.KASTOT 

50V. 

73 

continue 

510. 


GO TO (55,58*61) ,KOPT 

511. 

1000 

CONTINUE 

sw. 


RETURN , 

513. 


END 

514. 

♦ DATA 


515. 

INC 

45. 

516. 

€ P 

0.1 

517. 

HPO 

300. 50. 800. 

510. 

MASS 

too. 

519. 

0*G 

10. 

520. 

F10 

200. 

5?l « 

«srnp 


522. 

/• 
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